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We give results for the Upsilon spectrum from lattice QCD using an improved version of the 
NRQCD action for h quarks which includes radiative corrections to kinetic terms at 0{v^) in the 
velocity expansion. We also include for the first time the effect of up, down, strange and charm 
quarks in the sea using 'second generation' gluon field configurations from the MILC collaboration. 
Using the T 23 — IS splitting to determine the lattice spacing, we are able to obtain the IP — IS 
splitting to 1.4% and the 35' - IS splitting to 2.4%. Our improved result for M(T) - M{r]h) is 70(9) 
MeV and we predict M{T') — M(r]l) = 35(3) MeV. We also calculate tt, K and rjs correlators using 
the Highly Improved Staggered Quark action and perform a chiral and continuum extrapolation to 
give values for M^^ (0.6893(12) GeV) and (0.1819(5) GeV) that allow us to tune the strange 
quark mass as well as providing an independent and consistent determination of the lattice spacing. 
Combining the NRQCD and HISQ analyses gives mh/rris = 54.7(2.5) and a value for the heavy 
quark potential parameter of ri = 0.3209(26) fm. 



I. INTRODUCTION 

Lattice QCD calculations have developed rapidly both 
in accuracy and in scope in the last few years. This 
growth has built on the first demonstration that numer- 
ical simulations including li, d and s quarks in the sea 
with light enough u/d quarks give results in agreement 
with experiment for simple 'gold-plated' quantities across 
the full range of hadron physics [l] . Errors at the level of 
a few % make this highly non-trivial. A key element of 
those calculations was the determination of the T spec- 
trum because there are many gold-plated states below 
threshold for strong Zweig-allowed decay. In addition ra- 
dial and orbital excitation energies are very insensitive 
to quark masses (including that of the h itself) making 
them useful for determining the lattice spacing, a, with- 
out a complicated tuning process. A further incentive 
for lattice T studies is the importance of testing h quark 
physics from lattice QCD so that the same action can 
be used for results in B physics required, in conjunction 
with experiment, for the determination of elements of the 
Cabibbo-Kobayashi-Maskawa matrix. Here we give new 
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results for the T spectrum improving on those earlier re- 
sults in several ways to keep pace with improvements in 
other areas of lattice QCD. We have improved statistical 
errors, improved the NRQCD action and we are also now 
using 'second generation' gluon field configurations that 
include charm quarks in the sea. 

The h quarks in these first calculations that included 
the full effect of sea quarks [21 [3] were implemented us- 
ing lattice Nonrelativistic QCD (NRQCD) with an action 
accurate through in the velocity expansion for the h 
quark pfl . The coefficients of the terms were matched 
to full QCD at tree level, having removed the most signif- 
icant source of radiative corrections, that of tadpole dia- 
grams generated in lattice QCD from the form of the lat- 
tice gluon field, by the use of 'tadpole-improvement' [5]. 
The gluon field configurations used were generated by 
the MILC collaboration [6 using a Symanzik-improved 
gluon action in which radiative corrections at 0{asO?) 
were included except for radiative corrections from quark 
loops [7 {O{nfas0?) where is the number of sea quark 
flavors), which were omitted. Configurations at three dif- 
ferent values of the lattice spacing were available: 'super- 
coarse' (a ^ O.lSfm); 'coarse' (a ^ 0.12fm) and 'fine' (a ^ 
O.OQfm). u/d and s sea quarks were included using the 
improved staggered (asqtad) action [SHTO] which is nu- 
merically relatively fast. A range of u/d masses (taken 
to be the same) were used ranging down to a ratio with 
the s sea quark mass of around 0.2. The key mass split- 
tings in the bottomonium spectrum studied were those 
between the ground 6'-wave states and the first radially 
excited 6'-wave states, the 2S — IS splitting, and that 
between the first P-wave states and the ground 5'-wave 
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states, the IP — IS splitting. The statistical errors from 
the lattice calculation for these splittings were 1-2% (i.e. 
5-10 MeV) and systematic errors were estimated to be 
similar to this or smaller, depending on the lattice spac- 
ing. Within these errors, agreement with experiment was 
confirmed. 

More recently the T spectrum has been calculated us- 
ing the same NRQCD action on gluon configurations at 
a 'coarse' (a ~ 0.11 fm) and a 'fine' (a ~ 0.09 fm) lat- 
tice spacing generated by the RBC/UKQCD collabora- 
tion using the Iwasaki gluon action and 2+1 flavors of 
sea quarks implemented with the domain wall formal- 
ism [HI |T2] . Results in close agreement and with similar 
errors to those found on the MILC configurations are ob- 
tained, confirming the independence of the results from 
the sea quark formalism. 

The systematic errors in the calculation of the T 
25* — IS and IP — IS splittings were studied in some de- 
tail in [2j. Sources of error there were missing radiative 
corrections to the v"^ terms in the lattice NRQCD ac- 
tion (beyond tadpole-improvement), as well as radiative 
corrections to discretisation correction terms and from 
higher order (v^) missing relativistic corrections. In ad- 
dition systematic errors from the missing radiative cor- 
rections to the improvement terms in the gluon action 
were estimated. These errors were typically each of order 
1% in the IP — 16* splitting on the fine lattices and about 
half that for the 25* — IS* splitting because of some can- 
cellation between IS and 25* states. Errors were similar 
for the radiative and relativistic errors on coarser lattices 
but of course the discretisation errors were larger. 

Subsequent to this, we have made estimates of the ef- 
fect of missing c quarks in the sea [T^ |^ . These have 
negligible effect on mesons apart from bottomonium, 
where internal momenta can be large enough to gener- 
ate c quarks from the vacuum. We found the shift in the 
ground-state S wave masses might be of 0{5MeV) [M] 
(it is spin-independent) with approximately half the shift 
for 2S states because of a smaller 'wave-function at the 
origin' and no shift for IP states. This would give rise to 
systematic errors of 0.5% for the 2S — IS splitting and 
1% for the IP — IS splitting, similar to the systematic 
errors from other effects quoted above. 

The conclusion from these results is that the errors 
in bottomonium masses and radial and orbital mass 
splittings have been pinned down and tested from this 
NRQCD action at the level of 5-10 MeV. There is also a 
contribution to systematic errors at the same level com- 
ing from the gluon field configurations. The NRQCD 
systematic errors also feed in to the calculation of Bs 
and Be meson masses using NRQCD b quarks. The state- 
of-the-art calculation for the masses of these mesons has 
0(10 MeV) errors dominated by systematic errors from 
this NRQCD action [T1[15]. 

In the last five years, however, other lattice QCD cal- 
culations have become increasingly accurate. For exam- 
ple the mass of the Ds meson was recently calculated by 
HPQCD with combined statistical and systematic errors 



of 3 MeV and its decay constant calculated to 1% [13]. 
These errors are at the level where we must allow for 
missing electromagnetism from lattice QCD. 

There have been several contributions to this progress. 
Advances in computational speed have meant better sta- 
tistical errors from calculating many more meson correla- 
tors on larger samples of configurations. It has also been 
possible to generate lattices with smaller lattice spacing, 
so that the Ds calculation includes 'superfine' (a ^ 0.06 
fm) and 'ultrafine' (a ^ 0.045 fm) lattices [6 . Significant 
improvements have been made to relativistic quark ac- 
tions too. For example, the Ds meson mass calculation 
used the Highly Improved Staggered Quark (HISQ) ac- 
tion for both valence quarks. The HISQ action [16] has 
smaller discretisation errors than the asqtad action by 
about a factor of 3 and can be used for quarks as heavy 
as charm on lattices with a lattice spacing of O.lfm or 
smaller. This has revolutionised charm physics calcula- 
tions [17 in lattice QCD and is having an impact also 
on calculations for mesons containing a b quark through 
a combination of an extrapolations in the mass of the 
heavy HISQ quark acting as the '6' to the physical point 
for the real b quark, combined with extrapolations to the 
continuum (a 0) limit from results at many values 
of a [18]. The heavy HISQ calculations are computa- 
tionally much more expensive than those using NRQCD 
and this currently limits their utility. The results for Bs 
and Be meson masses have comparable errors to the ex- 
isting NRQCD results, but are dominated by statistical 
and a ^ extrapolation uncertainities. They then pro- 
vide a complementary way of testing b physics to that of 
NRQCD and it is clear that combining the strengths of 
both methods will be optimal in future. 

Meanwhile the MILC collaboration have moved on to 
the production of 'second generation' gluon field config- 
urations which have a number of improvements over the 
earlier ensembles [19]. They include a more highly im- 
proved gluon action [20^, HISQ quarks in the sea with 
the addition of c quarks as well as li, d and s and with 
lighter u and d masses than before. 

The availability of these configurations along with the 
incentives discussed above to improve errors in T and B 
physics using NRQCD b quarks has meant that we have 
begun a new programme of improved NRQCD calcula- 
tions. Here we present the first results, giving the radial 
and orbital splittings in the T spectrum, tuning the lat- 
tice b quark mass and determining the lattice spacing 
from the {2S — IS) splitting. As well as using the second 
generation gluon field configurations we have improved 
the NRQCD action by adding radiative corrections to 
the kinetic terms including discretisation errors. We 
also have improved statistics and improved methods for 
tuning the b quark mass. This has meant that we can 
test the effect of radiative corrections to the v"^ kinetic 
terms on the meson dispersion relation. Using both per- 
turbative and nonperturbative methods for determining 
the radiative corrections to spin-dependent terms we are 
able to improve the determination of the T hyperfine 
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splitting. 

A useful complementary method for determining the 
lattice spacing was developed in [21 . It uses the ficti- 
tious ss pseudoscalar particle known as the r]s • This par- 
ticle does not exist in the real world because of mixing 
with light quarks to form the 77 and rj' but on the lattice 
this can be prevented. The mass and decay constant of 
the can be determined accurately in a lattice QCD 
calculation using the Highly Improved Staggered Quark 
(HISQ) action and their physical values fixed from M^^, 
Mk, fn and /k from a simultaneous chiral and contin- 
uum extrapolation. Here we update the results of [21] for 
these 2-1-1+1 configurations and use these also to give a 
determination of the lattice spacing. 

The two different methods for determining the lattice 
spacing can be combined through the use of a third quan- 
tity, ri [22] , which can be derived accurately from deter- 
mination of the heavy quark potential [23 . ri/a values 
are provided for these configurations by the MILC col- 
laboration [24]. ri/a provides a good determinant of the 
relative lattice spacing between different sets of gluon 
configurations but its physical value must be determined 
from other quantities. From the separate determination 
of the lattice spacing from the two methods above we 
have two sets of results for ri in fm as a function of lat- 
tice spacing. From this we are able to test that the two 
methods give the same result in the continuum and chi- 
ral limits (which they do) and provide a physical value 
of ri that could be used, in the absence of either of the 
other methods, to determine the lattice spacing on other 
ensembles with 2+1+1 fiavors of sea quarks. 

We also combine results for tuned b quark masses in 
NRQCD and tuned s quark masses from HISQ along with 
one-loop renormalisation constants to give a value for 
miy/ms for comparison to other results obtained purely 
from the HISQ action. 

The layout of the paper is as follows. Section [H] dis- 
cusses the second-generation gluon field ensembles giving 
more details of the improvements present there. Sec- 
tion HI describes the improvements to the NRQCD cal- 
culations and results for the T spectrum. Section |IV| 
discusses the tt, rjs analysis on these same configu- 
rations and the additional information that provides to 
determine the lattice spacing. This is tied together via 
the determination of the heavy quark potential parame- 
ter, ri, in section[V|and m^/m^ in section VI Section VII 
provides our conclusions. 



II. SECOND GENERATION 2+1+1 GLUON 
FIELD ENSEMBLES 



The gauge configurations used in this calculation are 
listed in Table |l] [19 . These were generated by the MILC 
collaboration using a tadpole-improved Liischer-Weisz 
gauge action with coefficients corrected perturbatively 
through 0{as) including pieces proportional to tij, the 
number of quark ffavors in the sea [20 (see Appendix [A|). 



TABLE I: Details of the MILC gluon field ensembles used 
in this paper. /3 = 10/^^ is the SU{?>) gauge coupling and 
L/a and T /a are the number of lattice spacings in the space 
and time directions for each lattice, ami, arris and arric are 
the light (up and down taken to have the same mass), strange 
and charm sea quark masses in lattice units, ri/a is the static- 
quark potential parameter in lattice units determined by the 
MILC collaboration [191 EH- Note that this has not been 
'smoothed'. The ensembles 1 and 2 will be referred to in the 
text as "very coarse", 3 and 4 as "coarse" and 5 as "fine." 



Set 





ri/a 


ami 


ams 


amc 


L/a X T/a 


1 


5.80 


2.041(10) 


0.013 


0.065 


0.838 


16x48 


2 


5.80 


2.0621(45) 


0.0064 


0.064 


0.828 


24x48 


3 


6.00 


2.574(5) 


0.0102 


0.0509 


0.635 


24x64 


4 


6.00 


2.623(11) 


0.00507 


0.0507 


0.628 


32x64 


5 


6.30 


3.549(13) 


0.0074 


0.037 


0.440 


32x96 



The gauge action is then improved completely through 
0{asa'^), unlike the earlier asqtad configurations. Sea 
quarks are included with the HISQ action [16 which also 
has smaller discretisation errors compared to the asqtad 
action (see the discussion in section IV). The configura- 



tions include a sea charm quark in addition to up, down 
and strange. These configurations are then said to have 
2+1+1 fiavors in the sea, since the u and d quarks are 
taken to have the same mass, which is heavier than av- 
erage u/d mass in the real world, and the s and c masses 
are tuned as closely as possible to their correct values at 
that lattice spacing. The tuning of the sea s quark mass 
is much more accurately done - to better than 5% - than 
on the previous asqtad configurations. This means that 
the u/d quark mass (denoted mi here) can be more accu- 
rately calibrated in terms of the s quark mass for chiral 
extrapolations. Here we use a ratio of mi/rris as low as 
one tenth (see Table [l|) whereas in our previous work on 
the asqtad configurations our most chiral ensemble had a 
ratio of the mi^sea/'^s, physical of one quarter. This means 
that we have a much smaller chiral extrapolation to do 
to reach the physical u/d mass (where mi = ms/27 [6]) 
than before. 

The sea quarks are included with the standard method 
of incorporating the determinant of the quark matrix 
raised to the one quarter power for each fiavor, in or- 
der to implement the correct counting for sea staggered 
quarks. The algorithm used for including the sea quarks 
has now been improved by MILC to the exact RHMC 
algorithm [19 i.e. all errors in the time step for the up- 
dating algorithm have been removed. 

The configurations are separated by 5 trajectories in 
the time units of the updating algorithm for the very 
coarse and coarse ensembles and by 6 trajectories for the 
fine ensemble. In subsections IIIIBI and IIV Al we will 
study the autocorrelations in our meson correlators to 
show how independent the configurations are for different 
observables. 

The Ti/a values given in Table [l| are determined by 
the MILC collaboration after extraction of the potential 
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between two infinitely heavy (static) quarks at separation 
r/a in lattice units, ri/a is defined [22j as the point where 
the force F{r) derived from the derivative of the potential 
satisfies 



r'^F{r) = 1. 



(1) 



The values of ri/a for these ensembles have been chosen 
to match approximately those of the previous results in- 
cluding 2+1 flavors of asqtad quarks and can be used to 
determine the lattice spacing if the physical value for ri 
is known. Using the ri value determined previously on 
configurations with 2+1 flavors of sea quarks, this means 
that the lattice spacing values will be approximately 0.15 
fm, 0.12fm and 0.09fm. The physical spatial size of the 
lattices then exceeds 2.5 fm and is as high as 3.8 fm on the 
ensembles that correspond to mi/rris =0.1. In section [V| 
we will derive a physical value for ri based on the results 
from sections [TIT| and [TV| to calibrate more accurately the 
lattice spacing values for these configurations. 



III. THE UPSILON SPECTRUM 
A. The NRQCD action 

The spectrum of bottomonium mesons is extracted by 
computing appropriate correlators constructed from b- 
quark propagators on the gluon field ensembles listed in 
Table [l| We make use of NRQCD, an effective field theory 
that gives an expansion of the Dirac action in powers of 
the heavy quark velocity, v. This is discretised onto a 
space-time lattice as lattice NRQCD [H [25^ and is a 
good formalism to use for b quarks since they are known 
to be very nonrelativistic inside their bound states {v'^ ^ 
0.1). As used on the lattice NRQCD has the advantage 
that propagators can be generated using a simple time 
evolution equation rather than having to invert the Dirac 
matrix. The quark and antiquark fields are separated in 
this formalism as 2-component spinors. 

The NRQCD Hamiltonian we use is given by: 

aH = aHo + adH; 



adH 



2amb ' 



(A(2))2 



' 8{ami))^ 8(am5)^ x 
-C3 ,, i • (^V X E - E X 

-C4- Cr • B + C5 



V-E - E- V 



1 



-C6 



2amb 
(A(2))2 

16n(am5)2 * 



24am5 



(2) 



Here V is the symmetric lattice derivative and A^^^ and 
A^"^^ the lattice discretization of the continuum ^ • Df 
and Df respectively, arrib is the bare b quark mass. E 



and B are the chromoelectric and chromomagnetic fields 
calculated from an improved clover term [2 . The B and 
E are made anti-hermitian but not explicitly traceless, 
to match the perturbative calculations done using this 
action. 

In terms of the velocity expansion Hq is 0{v'^) and 5H 
is 0{v^)^ including discretisation corrections. Hq con- 
tains the bare quark mass parameter which is nonper- 
turbatively tuned to the correct value for the b quark as 
discussed below in subsection llll CI The terms in 6H have 
coefficients q whose values are fixed from matching lat- 
tice NRQCD to full QCD. This matching takes account 
of high momentum modes that differ between NRQCD 
and full QCD and so it can be done perturbatively, giv- 
ing the Ci the expansion 1 + c[^^as + 0{a1). In previous 
calculations [2 we used the tree level value of 1 for all 
the Q, after tadpole-improving the gluon fields. This 
means dividing all the gluon fields, Uj^{x) by a tadpole- 
parameter, uq^ before constructing covariant derivatives 
or E and B fields for the Hamiltonian above. The uq 
parameter corrects for tadpole diagrams that arise in a 
universal way from the way in which the lattice gluon 
field is constructed. For uq we took the mean trace of 
the gluon field in Landau gauge, uql- With tadpole- 
improvement in place we expect the radiative corrections 
to the Ci coefficients to be of normal size i.e. 0{1) [26]: 
without this they can be rather large. 

Here, on top of tadpole- improvement with uql, we use 
0{as) corrected coefficients for the kinetic terms, i.e. Ci, 
C5 and C6, so improving on the NRQCD action used previ- 
ously, and significantly reducing the systematic errors in 
the tuning of the b quark mass and in the determination 
of the radial and orbital mass splittings. The calculation 
of the c-^^ for z = 1, 5, 6 is discussed in Appendix [b| |27] . 
Table [ll| gives the values for ci, C5 and cq that we use on 
the very coarse, coarse and fine lattices as a result. As ex- 
pected, after tadpole-improvement, the coefficients c^^l g 
are not large and they are well-behaved as a function of 
the b quark mass. In subsection [III C| we test these coeffi- 
cients through a precision study of the dispersion relation 
for T and 7^5 mesons. 

The other coefficients in the NRQCD action are C2, 
C3 and C4. C3 and C4 multiply spin-dependent terms 
that give rise respectively to spin-orbit and spin-spin fine 
structure in the spectrum. Most of the splittings we will 
discuss here are 'spin- aver aged' to remove the effect of 
these terms and so we will generally set C3 and C4 to their 
tree level values of 1. However, in section [Til E 3| we will 
discuss the hyperfine splitting (M(T) — M{r]i))) and show 
results for both perturbatively improved and nonpertur- 
batively determined C4. The calculation of the appropri- 
ate is discussed in Appendix [b| and the non- 
perturbative determination of C4 and C3 in Appendix [C) 
The nonperturbative studies indicate that the value of 
C3 is very close to 1 for this NRQCD action. C2 multi- 
plies a spin-independent term, the Darwin term, which 
can affect spin-independent splittings such as radial and 
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TABLE II: The coefficients ci, C5 and cq used in the NRQCD 
Hamiltonian of equation [2] on the very coarse (sets 1 and 2), 
coarse (sets 3 and 4) and fine (set 5) ensembles. Other coeffi- 
cients had values 1 except for calculations in which we specif- 
ically changed their values to test the effect, as described in 
the text. 



Set 


Cl 


C5 


Cq 


very coarse 


1.36 


1.21 


1.36 


coarse 


1.31 


1.16 


1.31 


fine 


1.21 


1.12 


1.21 



orbital excitation energies. Because the Darwin term is 
field-dependent we do not expect it to have such a large 
effect as kinetic terms, and therefore do not expect ra- 
diative corrections to C2 to be as important as for ci, C5 
and Cq. However, in subsection [III C| we will investigate 
the effect of changing C2 so that we can estimate con- 
cretely the systematic error from not knowing its 0{as) 
correction. 

Given the NRQCD action above, the time evolution of 
the heavy quark propagator is given by: 



G(x,t + 1) = 1 



1 



a5H\ 
2n 



~2n 



uKx) 



1 



a6H 



with starting condition: 

G(x,0) = ,^(x)l. 



G(f,t)(3) 



(4) 



The smearing function ^(x) is used to improve the pro- 
jection onto a particular state in the spectrum. Includ- 
ing a variety of smearing functions is essential to obtain 
accurate results for the splittings between the low lying 
excited states. Full details of the smearing functions used 
will be given in subsection [III B| The 1 in equation |4] is 
the unit matrix in color and (2-component) spin space. 
The parameter n has no physical significance, but is in- 
cluded for improved numerical stability of high momen- 
tum modes that do not contribute to bound states [4 . 
In [2] it was demonstrated that radial and orbital mass 
splittings were the same within the statistical errors avail- 
able there for n = 2 and n = 4 on coarse lattices. The 
minimum value of n for stability increases as the h quark 
mass in lattice units falls on finer lattices. Rather than 
varying n as we change the quark mass, here we use n = 4 
throughout which is the value appropriate to the fine lat- 
tices. At zero spatial momentum the anti-quark propa- 
gator is the complex conjugate of the quark propagator 
for a source of the kind given in equation [4] 

Details of various parameters used in our calculation 



are listed in table III Tuning of the bare h quark mass will 
be discussed in subsection |III C[ The tadpole parameters 
uql were calculated by fixing a subset of each ensem- 
ble to lattice Landau gauge using a Fourier-accelerated 
steepest descents algorithm ^29^ to maximise the average 



TABLE III: Parameters used in the NRQCD action for our 
calculations that included a full 5x5 matrix of correlators. 
Other parameters have been used in subsidiary test calcula- 
tions as described in the text, arrih is the bare h quark mass 
and uql the Landau link tadpole- improvement factor used in 
the NRQCD action. The different number of digits given in 
the uql column reflect the precision with which it was deter- 
mined, ricfg gives the number of configurations used in each 
ensemble and nt is the number of starting time sources per 
configuration. Tp is the time length of each propagator in 
lattice units, asm is the parameter for the smearing function 
described in subsection IIII Bl 



Set 


arub 


UoL 


ricfg 


nt 


Tp 


dsm 


1 


3.42 


0.8195 


1021 


16 


40 


0.79 


2 


3.39 


0.82015 


1000 


16 


40 


0.80 


3 


2.66 


0.834 


1053 


16 


40 


1.0 


4 


2.62 


0.8349 


1000 


16 


40 


1.0 


5 


1.91 


0.8525 


874 


16 


48 


1.37 



trace link (Xl/i=i 4 ^ '^^^m(^))' which value, normalised, 
then becomes uql- The whole ensemble was then fixed to 
Coulomb gauge by using the same algorithm to maximise 
the spatial trace link (Xl^^i 3.2^ Tr[/^(x)) to allow us to 
use 'wave-function' smearing operators, with parameter 

Propagators were 



asm as described in subsection IIIB 
calculated from 16 time sources on each configuration to 
minimise statistical errors. Because in NRQCD we oper- 
ate a simple time evolution we can choose the time length 
of each propagator. This we take to be greater than or 
equal to half the time extent of the lattice as detailed in 
Table [ml 



B. Smearing functions and multiexponential fits 

Quark propagators are generated using three different 
smearing functions which we label as local, ground state 
and excited state. They are chosen to improve the pro- 
jection onto different radially excited states and previous 
experience has shown that 'hydrogen-like' wavefunctions 
work well [2]. 

(j^iir) = Sr,o 
(l)gs{r) = exp{-r/asm) 

(t)es{r) = (2as^ - r)exp(-r/a5^). (5) 

asm is the smearing radius and is chosen to be approx- 
imately the same in physical units for each ensemble. 
Values are given in Table [ITU Since a different smearing 
function can be applied separately to the quark and anti- 
quark we can make five different combinations as detailed 
in Table lEl 

A different smearing can also be applied at the source 
and the sink making correlator cominations labelled by 
e.g. lg,le,gG. The different smearing combinations allow 
the construction of up to a 5 x 5 matrix of correlators 
for the 5'-wave states that can be fit simultaneously. The 
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TABLE IV: Smearing combinations used for either the source 
or the sink in the construction of 5*- wave correlators. 



Name 


quark 
smearing 


ant i- quark 
smearing 


1 

g 

8 

G 
E 


(l)es 


4>i 

4>as 

4>es 



cross-correlators provide further useful information be- 
yond that in the diagonal terms that can be used in the 
fitting to extract the excited states more precisely. The 
correlators with quantum numbers of or ^S'o are dis- 
tinguished by the insertion of either a <j or a 1 in spin 
space at source and sink [30|. 

To make P-wave states we use only the 1 and g smear- 
ings above and apply a symmetric difference operator, 
A to the smeared source to give a P-'wavefunction'. 
This propagator is combined with that from a 5 func- 
tion source and a derivative applied at the sink to make 
a P-wave meson correlator. The complete set of com- 
binations of <j matrices with derivatives that are needed 
for the P-wave states is given in [30j . On the lattice the 
5-dimensional spin 2 representation is split into E and T2 
representations of the lattice rotational group and we fit 
these representations separately since differences in mass 
between them can arise from discretisation errors on the 
lattice. 

For the S'-wave states, statistical errors were improved 
further by using random wall sources in combination with 
the smearings discussed above. The delta function quark 
source is replaced with a (pseudo-) random colour vector 
ria{x) e U{1) at each spatial point of the initial time 
slice. When the meson correlator is constructed, the 
white noise property {r]a{x)r]\{y)) = 5ah^{x — y) ensures 
that the random noise cancels at all points except those 
where the initial spatial sites are the same. This can be 
combined with the smearing functions by distributing the 
random number associated with the centre of each smear- 
ing function along with the smearing function. Then once 
again the white noise property will mean that the resul- 
tant correlator averages over the initial time source the 
effect of having a smeared source at every point |31| . Pre- 
vious studies have found a significant improvement in the 
precision of the Upsilon ground state energy using ran- 
dom wall sources [ST. The improvement is less clear for 
excited states and therefore we did not use this technique 
for the P-wave states. 

Propagators were calculated from 16 time sources on 
each configuration but to avoid correlations between time 
sources, the correlators were binned over all sources on 
the same configuration. Autocorrelations between results 
on successive configurations in an ensemble were studied 



by calculating the autocorrelation function Cat 

(^i^i+AT) - (2:i)(Xi+AT) 



Cat — 



(6) 



Here Xi represents a correlator on a given ensemble, i. 
Xi-\-AT is the correlator on an ensemble separated by AT 
from i in the ordered ensemble i.e. AT = 1 corresponds 
to neighbouring configurations in the ensemble. The en- 
sembles have been generated taking into account the fact 
that autocorrelations increase on finer lattices. Thus 
neighbouring configurations are 5 trajectories apart for 
very coarse and coarse ensembles but 6 trajectories apart 
for the fine ensemble [19]. Cat is plotted against AT in 
Figure [1] for the case where x is an T correlator measured 
with a time separation on the lattice of approximately 
0.6 fm. This value was chosen to correspond to a point 
where correlators were dominated by the ground-state. 
The picture is qualitatively the same for different time 
separations, however. Cat drops to zero very rapidly, 
within the separation AT = 1. We therefore do not have 
to worry about autocorrelations between configurations 
but can treat them all as statistically independent. 

Bayesian fitting is used to extract the spectrum from 
the correlators 33 . The fit function 



5 ^s/c 



t) = a{nsc,k)a*{nsk,k)e~ 



■Ekt 



(7) 



is used, where aE^ is the energy of the {k — l)th ra- 
dial excitation in lattice units and a{nsc/skik) are the 
corresponding amplitudes labelled by the smearing used 
at the source and sink of the correlator, i.e. sc, sk G 
{l^g^e^G^E}. We fit the full range of t values for the 
correlator from 1 to T„ 



/S-wave fits in Table III 
number of terms, nexp, 



where Tp values are given for 



and Tp = 20 for P-waves. The 
in the fit is varied, however, and 
Bayesian model selection criteria are applied to deter- 
mine which fit is used. In practice, this means adding 
additional terms to the fit until the results and the er- 
rors stabilise. An example is given in Figure [2j 

The Bayesian approach allows the inclusion of prior 
data into the fitting procedure. The ^^st function is 
amended to 



Aaug 



2 I 2 

X + X 



prior 



(8) 



and the function Xaug minimised. By Bayes' theo- 
rem this corresponds to maximising the posterior prob- 
ability p (parameters I data) as opposed to a standard 
test which maximises only the likelihood function 
p(data|parameters). Xprior taken to be 



2 ^ {Pk -Pkf 

Xprior 



k 



(9) 



for each fit parameter p/^. This assumes that the prior 
probablility density function for each parameter is a 
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FIG. 1: Autocorrelation function Cat for T correlators made from different smearing combinations, from left to right: 11, gg 
and ee. Different symbols are given to different ensembles according to the key on the right in the ee plot (color online). The 
correlators are evaluated at lattice time separation t/a = 4 on very coarse lattices (sets 1 and 2), t/a = 5 on coarse lattices 
(sets 3 and 4) and t/a = 8 on fine lattices (set 5). This corresponds to a t value where the gg correlators have reached the 
ground- state plateau and the ee correlators have a short plateau corresponding approximately to the first excited state mass. 
AT gives the separation at which the autocorrelation is measured in units of numbers in the ordered ensemble list. 
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FIG. 2: Energies in lattice units of the low lying T states for 
the fine ensemble, set 5, from the full 5x5 IgeGE fit plotted 
against the number of exponentials, riexp, included in the fit. 



necessary to introduce a cutoff, Wcut^ on the lowest eigen- 
values of the correlation matrix in order to fit the data. 
A variation of this method is used in which, instead of 
setting eigenvalues below i^^cut'^max to zero, they are set 
to Wynax times Wcut- This is a less severe truncation of 
the correlation matrix and it improves the fits in some 
cases. Wcut was typically taken to be 10~^ for the 5x5 
matrix fits. 

In order to determine whether the inclusion of five dif- 
ferent smearing operators actually leads to improved re- 
sults, the energies of the low lying T states are plotted 
in Figure |3] for a variety of different matrix fits from the 
fine ensemble. The effect on the precision of the ground 
state is negligible but the full 5 x 5 fit has significantly 
smaller errors for the first two excited states. 

Because NRQCD is a nonrelativistic effective theory, 
there is an energy offset. Thus the energies obtained from 
correlators at zero momentum do not correspond to me- 
son masses. Energy differences do correspond to mass 
differences, however and so, for example, the mass differ- 
ence between the and the T (in lattice units) is given 
simply by aE2 — aEi from equation [7| To obtain absolute 
mass values requires the study of correlators for mesons 
at nonzero spatial momentum as discussed in Sec. |III C 



Gaussian with central value pk and width a^^. The fit 
parameters are: the amplitudes, which are taken to have 
a prior of 0.1 ±1.0; the ground state energies ln(£^o) which 
are estimated from an effective mass plot and given a suit- 
ably wide width; and the splittings ln(£^n+i — ^n) which 
prior information tells us should be of the order 500 MeV 
with a width of 250 MeV. Taking the fit parameters to 
be the logarithms of the energy splittings ensures that 
the ordering of the states is respected. 

Xaug is minimised using the singular value decomposi- 
tion (SVD) method. In the larger matrix fits, the cor- 
relation matrix can become ill-conditioned and it can be 



C. NRQCD systematics in tuning the b quark mass 

In this calculation the parameters of QCD that need to 
be determined are the b quark mass and Aqcd- In prac- 
tice this translates into the fact that we need to tune the 
b quark mass parameter in the lattice NRQCD Hamilto- 
nian until we obtain the correct value for one calibration 
hadron mass and we need to determine the lattice spac- 
ing from another calibration hadron mass. After that 
is done all other hadron masses are determined with no 
further tuning. The two calibration hadrons should be 
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1 Ig le Ige leE IgeGE 

FIG. 3: Comparison of the effect of different smearing combi- 
nations for extraction of the energies of the ground and first 
two radially excited T states. The energy in lattice units from 
the fine ensemble is shown for a 1 x 1 1 fit (plus), 2x2 fits Ig 
(star) and le (circle), 3x3 fits Ige (square) and leE (triangle), 
and the 5 x 5 fit containing all sources IgeGE (cross). 



chosen with rather different properties. The mass chosen 
to fix the b quark mass should ideally be very sensitive 
to that value; the mass chosen to determine the lattice 
spacing should be as independent of the b quark mass as 
possible to avoid a complicated iterative tuning process. 
To determine the lattice spacing we choose the radial 
excitation energy of the T, i.e. M{V) - M(T). This 
is known from experiment to be very insensitive to the 
heavy quark mass since it changes by only 4% between 
the b and the equivalent quantities for the c quark, which 
has a mass a factor of 4.5 smaller. The determination of 
the lattice spacing from this quantity will be discussed in 



section |III E| Here we focus on the tuning of the b quark 
mass and in particular on the effect of the improvements 
to the NRQCD action which we have implemented here 
for the first time. 



As discussed in section [TIIBI the fitted energy from a 
zero momentum hadron correlator made from NRQCD 
propagators is not the hadron's mass because there is 
an energy offset. Instead we must determine the 'kinetic 
mass' from the energy-momentum dispersion relation: 



,2 p2 



2aAE 



(10) 



where aAE is the energy difference between the meson 
with momentum Pa in lattice units and the meson at 



rest. Equation \iO\ assumes a fully relativistic dispersion 
relation, i.e. 



aE{P) = aE{0) + y/a^P^ ^ a^M^^. (11) 

Systematic errors will then be present in the kinetic mass 
for lattice NRQCD both because the action is only accu- 



rate to a specific order in the expansion in v^/c^ and from 
lattice discretisation errors. Here we study both of these 
effects. First it is worth briefly recapitulating a discus- 
sion from the literature (see, for example, [34 ) on how 
the kinetic mass is built up in a nonrelativistic approach 
as successive orders in v'^/c^ are added to the nonrela- 
tivistic expansion, because it provides a useful handle on 
systematic errors. 

By definition the mass of a meson is given by the sum 
of the masses of its constituent quarks plus the binding 
energy. The binding energy has contributions from the 
internal kinetic energy, i.e. the motion of the constituent 
quarks relative to the centre of mass, and from the inter- 
action energy. If we write the meson dispersion relation 
in the standard nonrelativistic expansion as: 



E{P) = Ml 



2Mo 



(12) 



then Ml is known as the static mass and M2 is the ki- 



netic mass, equal to M^^^ in equation 10 up to relativistic 
corrections. It should be possible to construct the cor- 
rect meson mass from both Mi and M2 i.e. the binding 
energy contribution needs to feed correctly into both of 
them. 

To see how this works in outline it is sufficient to study 
two free particles. The total energy of the two particle 
system is the sum of the masses, m^, plus the kinetic en- 
ergies, qf/2mi for each particle. In the center of mass 
frame (P = 0) this is simply mi + m2 plus the inter- 
nal kinetic energy. As is well-known, the internal kinetic 
energy can be written to leading nonrelativistic order as 
p^/2/i where p is the momentum of either particle in this 
frame and /i is the reduced mass (l//i = 1/mi + 1/7712). 
Thus Ml takes the expected form for this two particle 
system. To study M2 we must include the motion of the 
centre of mass and expand the sum of the two particle 
kinetic energies to O(P^). For M2 to have the correct 
form including the leading piece of the internal kinetic 
energy we need E{P) to take the form 



E{P) 



rriqi 



p2 



2/i 



(13) 



2{mqi +mq2) 



2/i(mqi +mg2) 



i.e. we need to locate a P^p^ term in the sum of the 
individual particle kinetic energies. This requires the in- 
dividual kinetic energies to be expanded beyond leading 
order in the nonrelativistic expansion to include terms 
at fourth order in the momentum. Thus M2 will have 
the correct form to leading order in the internal kinetic 
energy if the individual kinetic energy terms are correct 
through next-to-leading-order in momentum. In an in- 
teracting theory we also need the interaction terms to be 
correct through 0{v'^) to have the binding energy cor- 
rectly included in the kinetic mass. 

These issues are discussed in some detail in [34J for 
heavy quarks using the clover action since there are 
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FIG. 4: Kinetic mass values in lattice units obtained on the 
coarse ensemble, set 3, for the amb and Ci values given in 
Tables [ll] and |III| Kinetic mass values are given separately 
for the T and r/b and plotted against the square of the lattice 
momentum in units of 2na/L. The two results at x-axis value 
of 9 correspond to momenta with indices (3, 0, 0) and (2, 1,1). 
The higher one is (3,0,0). 



FIG. 5: Spin- averaged values for the kinetic mass in lattice 
units obtain ed o n the coarse ensemble, set 3, for arrib = 2.66 
(as in Table III). Results for the Ci values given in Table [u] 
are compared to the results for q = 1. The kinetic mass is 
plotted against the square of the lattice momentum in units 
of 2na/L. 



important differences in discretisation errors there be- 
tween choosing Mi or M2 as the appropriate meson mass 
against which to tune the quark mass. In NRQCD we 
must use M2 {M^-,^). The quark Hamiltonian given in 
equation |2] has no quark mass term, so to reconstruct 
the meson mass from Mi would require adding back in 
the zero of energy. This is perturbatively calculable but 
we wish the tune the quark mass fully nonperturbatively. 
M2 on the other hand acquires its quark mass pieces from 
the quark kinetic energy terms and so has no zero of en- 
ergy problem. As discussed above, M2 will also correctly 
include the internal kinetic energy if the relativistic 
corrections to the kinetic energy are included in the quark 
Hamiltonian, as they are in equation |2] Indeed we are 
now including the radiative corrections to the kinetic 
terms through adjustments to ci, C5 and ce, and we will 
show below the effect that this has. 

We can determine the kinetic mass very precisely by 
use of propagators made starting with a random wall 
source patterned by an exp(zp-x) factor to give the quark 
momentum [31]. We use only a 5{x) smearing function 
for these calculations so they are very fast, but we must 
evolve both a quark and an antiquark propagator be- 
cause the complex conjugate of a quark propagator of 
momentum p is an antiquark of momentum —p. Typi- 
cally we take quark and antiquark momenta to be equal 
so that the meson momentum, when they are combined, 
is P = 2p. 

We fit the meson correlator of momentum P simul- 
taneously with the meson correlator at rest so that the 
energy difference a/S.E between the ground state energies 
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FIG. 6: Spin- averaged values for the kinetic mass in lattice 
units obtained on the fine ensemble, set 5, for the arrib and Ci 
values given in Tables [ll] and [ill] compared to the results for 
Ci — 1. The kinetic mass is plotted against the square of the 
lattice momentum in units of 27ra/L. 



can be determined directly by the fit taking the correla- 
tions into account. In this way we obtain a/S.E values 
with errors typically in the 5th decimal place. To avoid 
cluttering the main body of the text, the detailed tables 
of values for T and 7^5 energies as a function of momentum 
and aMkin are collected in Appendix [D] Propagators were 
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calculated for the full number of configurations given for 
each ensemble in Table |III[ but in some cases we used 
fewer time sources per configuration than is given there. 

We can then plot out the kinetic mass for a range 
of meson momenta to study systematic effects in equa- 



tion [TO] which would show up as a disagreement between 
kinetic masses obtained from different momentum val- 
ues. Previous calculations saw no significant differences 
in kinetic mass values for momenta up to P^o? =9 with 
errors of around 1% |2j. This is equivalent to a test, as a 
function of momentum, of the constancy of the 'speed of 
light'. Here we are able to achieve errors down to 0.1%, 
depending on the momentum. Then systematic varia- 
tions of aMKin with momentum can be seen at the 0.5% 
level. 

aM^Kin values for T and 7^5 mesons on the coarse lat- 
tices, set 3, are plotted in Figure [4] and show several 
features. One is that there is a systematic difference be- 
tween the values of aM^,^ for on-axis (those in one lattice 
direction only) and off-axis momenta. This was hinted 
at in [2] but the errors were too large for it to be clear. 
The on-axis kinetic masses are higher, and this reflects a 
breaking of rotational invariance on the lattice which is 
a discretisation error. It is particularly obvious for the 
momenta with components along the spatial directions 
labelled by integers (3,0,0) and (2,2,1), both of which 
have P'^o? = 9(27ra/L)^. The difference is tiny but visi- 
ble. We will return to this point below. 

Another feature of Figure [4] is that the kinetic mass for 
the ?75 is above that of the T which is the opposite way 
round to the energy difference at zero momentum and 
to experiment. A similar but somewhat smaller effect is 
seen on the fine lattices. The discussion above on the 
way in which the meson kinetic mass is built up order by 
order in the nonrelativistic expansion shows how this has 
happened. It results from the fact that the a • B term 
that gives rise to the hyperfine splitting is only included 
at leading order in our NRQCD action, equation [2] Rel- 
ativistic corrections to this term would be needed for it 
to feed correctly into the kinetic mass, M2. The effect of 
the a ' B term splitting is correctly incorporated in the 
meson energy at zero momentum (Mi), however, and it 
is from differences in Mi for T and 775 that we determine 
the hyperfine splitting (see subsection III E 3 ) . This small 
but non-zero systematic error in M2 is simply removed 
by working instead with the spin-averaged kinetic mass 
of the T and 7^5: 



MKin(l^) 



(3Mk..(T)+Mk..(7?6)) 



(14) 



and using this to fix the h quark mass. 

The above arguments also allow insight into the effect 
of radiative corrections to the kinetic terms in the 
NRQCD Hamiltonian that we include here for the first 
time. Changing the coefficient of the p^/8m^ term, ci, 
from 1 to 1 + 0{as) will modify the amount of the inter- 
nal kinetic energy that is incorporated into the meson ki- 
netic mass, effectively correcting for an 0{as) mismatch 



between this contribution to Mi and M2 from binding en- 
ergy effects. The effect of this radiative correction is seen 
clearly in Figure [5] where we compare the spin-averaged 
kinetic mass with all q set to 1 to that from having the 
radiatively improved coefficients given in Table [TTj The 
difference would be expected to be 0{as x B) where 
B is the binding energy of (9(500 MeV). This could in 
principle be as large as 150-200 MeV. From Figure [5] we 
see that the effect is somewhat smaller than this on the 
coarse ensemble set 3 - a shift of kinetic mass of 0.05 in 
lattice units corresponds to around 80 MeV on these lat- 
tices. The shift is clearly visible, however. The radiative 
correction acts to increase the kinetic mass for a given 
bare h quark mass. This is because ci > 1 and the bind- 
ing energy is positive. Thus the correctly tuned quark 
mass will be lower (by the same percentage shift as that 
for the kinetic mass) when radiative corrections are in- 
cluded. A similar shift is observed on the fine lattices as 
shown in Figure [6j 

Remaining systematic errors from higher order radia- 
tive corrections to terms in the NRQCD action will 
be suppressed by a further power of as beyond the shift 
seen here. We therefore expect the remaining error in the 
kinetic mass from this source to be 0(0.3%). System- 
atic errors from missing higher order, ^ terms at tree 
level in the NRQCD action are a factor of v'^ ^ or 10%, 
smaller than the size of the effect of terms, and there- 
fore of similar size to missing a^v^ terms. They will also 
have the effect of correcting for momentum-dependence 
in MKin- From Figure [5] we can see that there is a sign 
of an upward drift of M^^^ with momentum but the ef- 
fect is smaller than the shift of Mxin with the radiative 
correction to the q coefficients. 

We now return to the issue of discretisation errors in 
the kinetic mass. These arise from the replacement of 
time and space derivatives in the NRQCD action with fi- 
nite differences on the lattice. The terms with coefficients 
C5 and cq contain o?v^ and av^ correction terms to re- 
move these errors. With the inclusion of radiative correc- 
tions to C5 and C6, the remaining errors are at 0{alo?v^) 
in this calculation. The term with coefficient C5, i.e. the 
term proportional to A^^^ is of interest because this is ro- 
tationally non- invariant. The signal for a lack of contin- 
uum rotational invariance in our results is a disagreement 
between the kinetic mass for on-axis momenta, that typ- 
ically have a high value for P/, and off- axis momenta. 
This was seen in Figure [4] for the coarse lattices. Less 
variation is evident on the fine lattices (Figure [6| , as ex- 
pected for a discretisation effect. 

To make clearer the way in which the rotationally non- 
invariant discretisation errors depend on the lattice spac- 
ing Figure [7| plots the energy difference in physical units 
between mesons with momentum (3,0,0) and (2,2,1) as a 
function of o? using results from all three values of the 
lattice spacing. P'^o? = 9(27ra/I/)^ corresponds to ap- 
proximately the same physical momentum at all three 
lattice spacing values, so the results should be a good 
test of how rotational invariance is restored as a ^ 0. In 
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FIG. 7: The energy difference in MeV between mesons with 
momentum (3,0,0) and (2,2,1) in units of 2na/L on the lattice 
plotted against the square of the lattice spacing in fm. Results 
are shown for the case ci,5,6 = 1 as well as for ci,5,6 taking 
their as-improved values. An example fit curve with and 
dependence is shown through the T data for ci,5,6 ols- 
improved. 



fact the energy difference is tiny on all except the very 
coarse lattices, where it reaches 1 MeV. The case in which 
the ci^5^6 coefficients are set to their tree level values of 
1 is plotted as well as the case with the ci^s^e coefficients 
taking the radiatively improved values that we have used 
for the rest of our calculation here. The radiatively im- 
proved values give very slightly smaller energy splittings, 
since they have improved the contribution to this er- 
ror by one order in as to 0^0? v^. The energy difference 
between mesons with momentum (3,0,0) and (2,2,1) also 
has contributions at 0{a^v^)^ however, and both the ef- 
fect of radiative improvement and the shape of the curve 
in Figure [t] tend to imply that these terms dominate 
over any remaining terms. 

Rotationally invariant discretisation errors would give 
rise to a kinetic mass that varied with . This is the 
same effect as that of relativistic errors, because the cor- 
recting operators are the same. Discretisation errors re- 
quire an a-dependent coefficient to correct them. How- 
ever, as discussed above under relativistic corrections, 
there is no sign in our results of such errors to better 
than 0.5%. 

The conclusion from this subsection is that, to min- 
imise systematic errors, we should tune the h quark mass 
by calculating the spin- averaged kinetic mass MKin(l^) 
and matching that to experiment. We do this from the 
comparison of meson energies at zero momentum and 
the 'maximally off- axis' momentum (1,1,1) to minimise 
discretisation errors. Table [V| gives results for this ki- 
netic mass on all ensembles for the given values of the h 
quark mass and coefficients, q. To convert these results 
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FIG. 8: Comparison of values obtained for the kinetic mass 
from a variety of different parameter values on coarse set 3. 



to physical units we need a value for the lattice spacing 
to be determined in subsection |III E| Table [V| gives sta- 
tistical/fitting errors on the values. As discussed above, 
remaining systematic errors from missing radiative, rel- 
ativistic and discretisation errors amount to a total of 
0.5%. We are able to pin down the size of these system- 
atic errors by using the improved methods described here 
to study the dispersion at this level of detail. 

Figure [8] compares the results for the spin-averaged ki- 
netic mass on the coarse ensemble, set 3 for a variety of 
different choices for the coefficients in the NRQCD action 
to show the size of variations in the kinetic mass. The 
figure shows that we can see the difference between tak- 
ing tree-level values for ci^s^e and radiatively improved 
values. Changing C2 (the coefficient of the Darwin term) 
has very little effect. The effect of changing C4 (the co- 
efficient of the cr • B term which should be spin-averaged 
away at leading order in this kinetic mass) is also not 
large. Another check of this is given in Table [V| on set 1. 

The experimental result for the T mass is 9.4603(3) 
GeV and that of the 7^5, 9.391(3) GeV, [35 giving a spin- 
average of 9.443(1) GeV. The real world includes effects 
that are missing from our lattice calculation, however, 
and so we must correct for this. Electromagnet ism affects 
the T and 775 approximately equally and, from a potential 
model we estimate that it reduces their masses by 1.6 
MeV [14]. In addition the 7^5 can annihilate to gluons 
and we estimate that this effect also reduces its mass by 
2.4 MeV, taking the same value as that estimated for the 
Tjc [16 . The 'experimental' mass that we should compare 
our results to is then increased from above to 9.445(2) 
GeV where we allow for a 100% error in our estimate of 
the shifts in the masses 60 . 
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as 


1.25 1 5.748(8) 


5.823(4) 


5.766(7) 


3 


2.66 


as 


1 1.25 5.767(10) 


5.889(4) 


5.798(8) 


4 


Z.DZ 


as 


1 1 5.706(9) 


5.761(4) 


c Ti ^\^'^7\ 

5.719(7) 


4 


2.66 


as 


1 1 5.778(11) 


5.833(5) 


5.792(10) 


5 


1.91 


1 


1 1 4.230(13) 


4.252(6) 


4.236(11) 


5 


1.91 


as 


1 1 4.256(14) 


4.287(6) 


4.264(11) 


5 


2.0 


as 


1 1 4.431(11) 


4.466(5) 


4.439(10) 



TABLE V: Summary of the kinetic masses obtained on different ensembles for a variety of parameter values. We use the energy 
difference between lattice momentum zero and momentum ap = (1, 1, 1) in units of 27ra/L. The column ci,5,6 denotes whether 
the 0{as) improved coefficients were used in the action and the columns C2,C4 indicate additional values of those coefficients 
that were run on coarse set 3 and very coarse set 1 to estimate systematic errors. 



D. NRQCD systematics in radial and orbital 
splittings 

Here we discuss the remaining sources of systematic er- 
ror in our calculation of the radial and orbital excitation 
energies. These systematic errors will feed subsequently 
into the determination of the lattice spacing from the T 
2S -IS splitting. 

Radial and orbital excitation energies arise at leading 
order from the time derivative and Hq in the NRQCD 
action (equation |2|. The relativistic corrections at 
in SH thus provide relative 0{v'^) ~ 10% corrections to 
these splittings. Missing radiative corrections to the v"^ 
terms dominated the errors in earlier calculations [2j [11] , 

since agv'^ ~ 2 3% is larger than ^ 1% from missing 

higher order relativistic corrections. We now include for 
the first time the radiative corrections to most of the 
terms in SH. The remaining errors are then largely at 
relative 0(a^'i;^), i.e. less than 1%. 

Table |VI| lists the remaining systematic errors from 
spin-independent terms in the 25* — IS and IP — IS 
splittings in more detail following [2 . The errors were 
determined using a potential model to make estimates 
of the energy shifts in each of the IS^ 2S and IP 
states. For example, radiative corrections at 0{al) to 
the p^/(8m^) term in the NRQCD action give shifts of 
size al < > /4m^ where < > is the expectation 
value of in that state. 

The effects of the Darwin term term appear at 0{asV^) 
since we have not included a radiative correction to C2. 
However, since this term vanishes in the free theory it 
is already suppressed by an additional power of a^. Its 
effects are proportional to the square of the wavefunc- 
tion at the origin so it does not affect P-wave states. A 
very similar term arises from missing spin-independent 
4-quark opertors. The spin-dependent 4-quark operators 
are discussed in Appendix [B] along with the coefficients 
they have in order to match NRQCD to QCD. The spin- 



independent ones arise from the same diagrams and the 
calculation of their coefficients is in progress. Here we 
take an error from missing these 4-quark operators which 
is of the same size as the error from radiative corrections 
to the Darwin term. 

Note that errors cancel to a significant extent between 
the 25* and IS* states because of their similarities [2j. This 
is the reason for focussing on the 25* — IS* splitting to 
determine the lattice spacing, because it has the smallest 
systematic error. 



We see from Table VI that the largest remaining sys- 
tematic error is now that from missing terms. The 
key kinetic term at that would appear in a higher or- 
der NRQCD action is -(A(2))^/(16(am5)^) at tree level. 
This term is proportional to and so, if it domi- 

nates the errors, they will have the same sign at ev- 
ery value of the lattice spacing. Including this term 
would act in the direction of reducing both the 2S — IS 
and IP — IS splitting but the IP — IS splitting would 
be reduced the most. 



Table VII similarly quantifies remaining systematic er- 
rors from missing al radiative corrections to the discreti- 
sation correction terms with coefficients C5 and cq. These 
are significantly reduced over our earlier calculations ^ 
now that the radiative corrections are included. In 
addition the gluon action is now improved completely 
through 0{alo?) [20 and this means that the discreti- 
sation errors coming from the gluon action are similarly 
reduced. 

We can estimate the size of errors from the analysis 
in subsection |III C| where we study discretisation errors in 
the kinetic mass. The energy difference between mesons 
of momenta (3,0,0) and (2,2,1) in units of 27ra/L can 
be taken as a measure of at least the rotationally non- 
invariant errors, as discussed there. The energy dif- 
ference (Figure [7|) is barely visible except on the very 
coarse lattices where it amounts to 1 MeV, or 0.2% of 
the 25* — 15* splitting. This is much less than the esti- 
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Correction 




relativistic 


radiative 
kinetic 


radiative 
Darwin 


4- quark 
spin-independent 


Total 

relativistic 
+ radiative 




Form 














Est. 


%age in 25' - 


IS 














very coarse 




0.5 


0.2 


0.4 


0.4 


0.8 




coarse 




0.5 


0.15 


0.3 


0.3 


0.7 




nne 




U.O 


0.1 


O.Z 


O.Z 


O.D 


Est. 


%age in IP — 


IS 














very coarse 




1.0 


0.7 


0.9 


0.9 


1.8 




coarse 




1.0 


0.5 


0.7 


0.7 


1.5 




fine 




1.0 


0.3 


0.4 


0.4 


1.2 



TABLE VI: An estimate of systematic errors in the 2S — IS and IP — IS splittings in the T in our lattice QCD calculation 
arising from missing higher order relativistic and radiative corrections to the NRQCD action that we use (equation [21 . 





Correction 




discretisation in 


discretisation in 


discretisation in 


Total 








NRQCD action (i) 


NRQCD action (ii) 


gluon action 


discretisation 




Form 




aiaSp"^ / Sn{mby 


aia^dpi /12mb 


47ra^a^^(0)Vl5 




Est. 


%age in 2S - 


IS 












very coarse 




0.2 


0.4 


0.3 


0.5 




coarse 




0.1 


0.2 


0.15 


0.3 




fine 




0.05 


0.06 


0.05 


0.1 


Est. 


%age in IP — 


IS 












very coarse 




0.7 


2.0 


1.0 


2.3 




coarse 




0.4 


1.0 


0.5 


1.2 




fine 




0.2 


0.3 


0.1 


0.4 



TABLE VII: An estimate of systematic errors in the 2S — IS and IP — IS splittings in the T in our lattice QCD calculation 
arising from discretisation errors in the NRQCD and gluon actions. 



mate of remaining errors in that case so we do not 
include it in Table [VTll 



E. Results 

1. Radial and orbital excitation energies 

Our main results for the fitted energies for the ground- 
state and first two radial excitations of the T and 7^5 



are given in Table VIII The values come from multi- 
exponential fits to a 5 X 5 matrix of correlators for each 



meson as described in section IIIB We take 9 exponen- 
tials on sets 1, 2 and 3; 11 exponentials on set 4 and 12 
on set 5. We also give the fitted ground-state energy for 
the h})(lP) state on sets 3 and 5 from a 5 exponential fit 
to 2 X 2 matrix of correlators. The h quark masses and 
coefficients, used in the NRQCD action are those of 
Tables [ll] and III Errors are very small on the ground- 
state ^-wave masses but increase rapidly with the radial 
excitation number. The table also includes energy split- 
tings in lattice units for radial and orbital excitations. 

As explained earlier we can use the radial excitation 
energy, M{T') — M(T), to fix the lattice spacing, by 
setting 



a-^(GeV) = 



0.5630(9) 



aE(2^Si)-aE(l^Si) 



(15) 



0.5630(4) GeV is the experimental mass difference and 
we have increased the error to allow for a possible rela- 
tive shift in the two masses as a result of the electromag- 
netic attraction between quark and antiquark missing in 
our calculation. As discussed earlier, a potential model 
estimate would give a shift of 1.6 MeV to the T from 
the electrostatic attraction between quark and antiquark, 
and somewhat less for the since typical separations be- 
tween quark and antiquark are larger. We do not shift 
the result but allow for an error of 0.8 MeV. 

As long as we deal with spin-averaged splittings we 
do not have to consider errors in spin-dependent terms. 
However, for the 28 — IS splitting the match to experi- 
ment cannot be spin-averaged since no experimental in- 
formation is available for the 775(26'). In that case we have 
to consider sources of systematic error in the hyperfine 
splitting that will induce errors in the T and energies. 
This will discussed further in subsection IIII E 3l 

The main source of error is from missing radiative cor- 
rections when we take the coefficient of the a • B term, 
C4, to be 1. In section |IIIE3| we compare results for 



C4 = 1 to those from C4 corrected perturbatively through 
0{as) and nonperturbatively, to give the correct l^P 
fine structure. Both methods for correcting C4 give val- 
ues above 1 and increase the lattice result for the hy- 
perfine splitting (which is proportional to c\ at leading 
order). Thus with C4 = 1 the T energy is too low. Since 
M(T) = M(lS^) + (M(T)-M(7^6))/4, the shift from C4 in 
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aE{l'So) 
aE{2^So) 
aEiS^So) 
aEil^'Si) 
aEi^'Si) 
aEiS^'Si) 



0.25080(5) 

0.6898(16) 

0.975(14) 

0.28532(6) 

0.7078(14) 

0.988(16) 



0.25361(3) 

0.6909(8) 

0.940(22) 

0.28809(3) 

0.7074(8) 

0.975(8) 



0.26096(3) 

0.6235(8) 

0.849(9) 

0.29245(3) 

0.6416(7) 

0.855(11) 



0.26524(2) 

0.6246(6) 

0.854(4) 

0.29681(2) 

0.6393(9) 

0.867(10) 



0.25851(2) 

0.5248(7) 

0.677(11) 

0.28405(2) 

0.5370(9) 

0.693(10) 



aE{2S - IS) 
aEiSS - IS) 
aEl2^So - l^So) 
aEiS^So - l^So) 
aEl2^Si - l^Si) 
aEiS^Si - l^Si) 



0.4266(11) 

0.708(12) 

0.4390(16) 

0.724(14) 

0.4225(14) 

0.703(16) 



0.4238(7) 

0.687(8) 

0.4373(8) 

0.687(22) 

0.4193(8) 

0.687(8) 



0.3525(6) 

0.569(9) 

0.3626(8) 

0.588(9) 

0.3492(7) 

0.563(11) 



0.3467(7) 

0.575(8) 

0.3594(6) 

0.588(4) 

0.3425(9) 

0.570(10) 



0.2563(7) 

0.411(8) 

0.2663(7) 

0.418(11) 

0.2530(9) 

0.409(10) 



Rs 
aA 




1.664(38) 1.638(19) 1.611(32) 
0.00190(1) 0.00190(1) 0.00151(1) 


1.665(31) 1.617(40) 
0.00151(1) 0.00091(1) 


aE{l'Pi] 

aE(l^Pi 

Rp 


-IS) 


0.5654(23) 
0.2809(22) 
0.808(7) 


0.4833(10) 
0.2056(10) 
0.816(5) 


aEil-'Si 
aEl2^Si 
Rh 


-I'So) 
-2' So) 


0.03452(8) 0.03448(4) 0.03149(4) 0.03157(3) 0.02554(3) 
0.0180(21) 0.0165(11) 0.0181(11) 0.0147(10) 0.0122(11) 
0.521(62) 0.479(33) 0.575(35) 0.465(34) 0.478(45) 



TABLE VIII: Radial, orbital and 5'-wave fine structure splittings in lattice units for sets 
and coefficients given in Tables [TT| and III cs 
text. 



C4 = 1.0. Errors are statistical/fitting only. 



1 to 5 for the NRQCD parameters 
Rs, Rp and Rh are defined in the 



C3 

C4 



3 

1.0 

1.25 



4 

1.0 
1.25 



5 

1.0 
1.10 



aE{l'So) 
aE{2^So) 
aE{3^So) 
aEil^'Si) 
aE(2''Si) 
aEiS^'Si) 



0.20943(3) 

0.5796(6) 

0.788(12) 

0.25628(4) 

0.6022(7) 

0.827(7) 



0.21289(2) 

0.5777(7) 

0.802(6) 

0.25978(2) 

0.5999(7) 

0.821(5) 



0.23204(2) 

0.5021(12) 

0.660(12) 

0.26206(3) 

0.5170(18) 

0.663(24) 



aE{2S - IS) 
aElsS - IS) 
aEl2^So - l^So) 
aEiS^So - l^So) 
aEl2^Si 
aEls^Si 
Rs 



I'S^) 
I'S^) 



0.3520(6) 

0.573(6) 

0.3702(6) 

0.579(12) 

0.3460(7) 

0.571(7) 

1.651(20) 



0.3463(6) 

0.568(4) 

0.3648(7) 

0.589(6) 

0.3401(7) 

0.561(5) 

1.650(15) 



0.2584(13) 

0.405(18) 

0.2701(12) 

0.428(12) 

0.2549(18) 

0.401(24) 

1.573(95) 



aE{l'Pi] 




0.5247(22) 0.5253(20) - 


aE{l^Pi 


-IS) 


0.2801(22) 0.2773(20) - 


Rp 




0.810(7) 0.815(6) - 


aE(l^Si 


- I'So) 


0.04684(5) 0.04689(3) 0.03003(4) 


aEl2^Si 


-2' So) 


0.0226(9) 0.0222(10) 0.0149(22) 


Rh 




0.482(19) 0.473(21) 0.496(73) 



TABLE IX: Radial, orbital and 5'-wave fine structure splittings in lattice units for sets 3, 4 and 5 with NRQCD coefficients 
and parameters as in Tables [ll] and \nj\ In addition C4 is nonperturbatively tuned taking values from Table pClHl C3 = 1 for all 
results. Errors are statistical/fitting only. Reduced statistics of 400 configurations were used for the 5'-wave states from set 5. 
Rs, Rp and Rh are defined in the text. 



the T mass is one quarter of the change in the hyperfine 
splitting. In section |IIIE3| we also determine the ratio 
of the 2S hyperfine splitting to that of the 15* hyperfine 
splitting and find a result close to 0.5, independent of C4. 
Thus the shift from a change in C4 to the 2S — IS split- 
ting is one eighth of change in the IS hyperfine splitting. 



An increase in C4 above 1 reduces the 2S — IS splitting. 
From Table |XlV| we can compare results for the IS hy- 
perfine splitting for C4 = 1 to the value obtained for C4 
improved through 0{as) for sets 1, 3 and 5. Dividing 
by 8 then gives shifts in lattice units that can be applied 
to correct the 2S — IS splitting on very coarse, coarse 
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Set ar (fm) 



an 



(fm) 



1 0.1474(5)(14)(2) 0.1546(10)(5) 

2 0.1463(3)(14)(2) 0.1526(6)(5) 

3 0.1219(2)(9)(2) 0.1234(7)(4) 

4 0.1195(3)(9)(2) 0.1218(5)(4) 

5 0.0884(3)(5)(1) 0.0899(6)(3) 



ar^/a (fm) 

0.1569(8)(13) 

0.1553(3)(13) 

0.1244(2)(10) 

0.1221(5)(10) 

0.0902(3)(7) 



TABLE X: Lattice spacing values in fm determined from sev- 
eral methods. The first column gives results from the T 
25'— IS splitting. The first error is from statistics/fitting, the 
second from remaining sys tema tic errors from the NRQCD 
action (from Tables [Vl| and VII) and the third is a correlated 



0.2% error from experiment and electromagnetic corrections. 
The second column gives lattice spacing values from the decay 
constant of the ris meson as described in section [lV| The first 
error is from statistics/fitting and the second is a correlated 
0.3% error from the uncertainty in the physical value of fr^^ as 
discussed in section [lV| The third column gives lattice spac- 
ing values determined from ri/a values in Table [l] The first 
error is from statistics/fitting and the second is a correlated 
0.8% error from the uncertainty in the physical value of ri as 
discussed in section [Vl 



and fine lattices. These shifts are denoted by aA in Ta- 
ble |VIII[ and are to be subtracted from the 2S' — IS* split- 
ting to give the corrected lattice result. It can be seen 
that aA is not much larger than the statistical errors on 
the 25* — IS splitting. The statistical error in aA is neg- 
ligible, but there is a systematic error which is taken as 
0.5 X a A. This accounts for the errors in the hyperfine 
splitting from 4-quark operators, higher order radiative 
corrections to C4 and relativistic corrections to the a • B 
term. This error is then included in the systematic error 
for the corrected 2S — IS splitting. 

Note that we do not expect the spin-orbit term with co- 
efficient C3 to have significant effect on the 5'-wave states. 
In any case our nonperturbative determination of C3 dis- 
cussed in Appendix [C] gives a result consistent with the 
value of 1.0 that we are using. Possible errors from ra- 
diative corrections to C2 are included in our systematic 
error budget for NRQCD (Table |Vl|. 

Table [X] gives the values of the lattice spacing in fm 
obtained from the 25* — 15* splitting on each ensemble, 
along with their associated statistical /fitting error and 
systematic error. The systematic errors are combined in 



quadrature from Tables VI and VII and from aA in Ta- 



ble VIII The systematic errors are dominated by those 
from missing higher order relativistic corrections to the 
NRQCD action and these will be correlated to some ex- 
tent between ensembles. There is an additional overall 
systematic error of 0.2% coming from the experimental 
value for the splitting and electromagnetic effects missing 
from our calculation. 

In Table [TX] we give results for cases where C4 is set to 
its nonperturbatively tuned value on sets 3, 4 and a test 
value of 1.10 on set 5 (the nonperturbatively tuned value 
is in fact 1.18, see Appendix [C|) . Changing C4 shifts the 
fitted energies of all the states but this is simply because 
the zero of energy has changed. As expected, changing C4 
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FIG. 9: Results for the ratio of the 35'- 15* and IP- IS' split- 
tings to the 2S—1S in the T system plotted against the square 
of the lattice spacing determined from the 2S — IS splitting. 
The grey shaded bands give the physical result obtained from 
a fit to the data as described in the text. The black open 
circles slightly offset from a — are from experiment 35 . 



has very little effect on splittings between spin-averaged 
S wave states or between the ^Pi mass and the spin- 
averaged 15* state. 

An important test of the results is whether, using these 
values for the lattice spacing, we get results in agreement 
with experiment for other mass differences i.e. whether 
ratios of splittings are correct. In our previous work on 
2+1 flavor gluon configurations [2 agreement with exper- 
iment was found within 3% statistical/systematic errors. 
Here we have substantially improved errors, including im- 
proved statistical errors, so we can improve on our earlier 
analysis. 

Table [VIII gives values for the ratios of the T 3S — IS 
and l^Pi — 15* splittings to the T 2S — IS splitting from 
our results for C4 = 1 . Table [IX] gives the same ratios for 
the case where C4 takes its nonperturbatively tuned value. 
In forming the ratio Rp = (l^Pi - IS)/ {2^ Si - l^Si) for 
the case C4 = 1 we correct the denominator for C4 er- 



rors using the a A values in Table VIII The numerator 
should not be sensitive to C4 because the S'-state energies 
have been spin- averaged and the ^Pi state is unaffected 
by C4, as discussed in Appendix [C) Both the T 36* — 15* 
and 2S — IS splittings will have some sensitivity to C4. 
However, any shifts will cancel between the two splittings 
up to an amount equal to one quarter of the difference 
in the 35* and 2S' hyperfine splittings. This is negligi- 
ble compared to statistical errors in this ra tio. T he ratio 
Rs = {3^ Si - 1^ Si)/ {2^ Si - l^Si) in Tablef 



VIII 



^is there- 
fore not corrected for C4. 

The ratio Rs from Table VIII[ where we have results 
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for all five sets, is plotted against the square of the lattice 
spacing in Figure [9] We see very little dependence on 
lattice spacing and on the u/d sea quark mass (the s and c 
sea quark masses are already well-tuned to their physical 
values, see section IV). Figure [9[also shows results for i?p, 
combining results from Tables VIII and IX since we have 
results for only 2 ensembles in each table. The results 
change very little between the ensembles, however, as 
Figure [9] shows. 

To derive a physical value for each ratio we can use the 
results to fit for the dependence on these two quantities 
and then determine the result at physical sea quark mass 
and at a = to compare to experiment. For the sea u/d 
quark mass dependence a simple polynomial in mi/rris 
(values in Table |l| suffices because the mi values are al- 
ready very close to the physical point with mi/rris = 0.1 
and 0.2. The dependence on the lattice spacing is more 
complicated because in NRQCD we must allow for un- 
physical a-dependence coming from am^-dependent ra- 
diative corrections to discretisation errors. This am^- 
dependence is mild when ami) is sufficiently large, as here, 
and this is seen explicitly in the radiative corrections that 
are included in our calculation for C5 and cq (Table |llj). 

We therefore fit each ratio, i?, to the following func- 
tional form: 



These values are plotted along with the lattice results in 
Figure |9] 

The complete error budget for the two ratios is given in 



R — Rphys [1 



(16) 



+ 2bi5xi{l ^ CliaAf)] . 

Here Sxi is {ami/ am s) — (^^/^s)phys foi" each ensemble. 
{mi/ms)Y>h.ys is taken from lattice QCD as 27.2(3) [6]. 
The strange sea quark mass is tuned to better than 3% 
with the lattice spacing taken from the T 2S — IS split- 
ting so we can ignore any effects from this mistuning 
since sea quark mass effects are so small. 5xm allows for 
variation in the value of am^ over the range we are using 
and therefore a change in the NRQCD radiative correc- 
tions to discretisation errors. We choose 5xm to vary 
from -0.5 to +0.5 over our full range of masses by setting 
^^m = {ami) — 2.65)/1.5. A sets the scale for physical 
a-dependence. We take it to be 500 MeV. 

The fit prior on Rphys is taken to be 0.8 ± 0.1 for Rp 
and 1.6 ±0.2 for Rs. Since tree-level errors have been 
removed in this calculation we take the prior on terms 
to be 0.0 ±0.3; we take 0.0 ±1.0 for higher order terms in 
a. For bi we take 0.0 ±0.015 allowing for a 3% shift if the 
u/d quarks were as heavy as strange. Previous results [2] 
saw a 10% shift in results in the quenched approximation. 

Good fits using the form in equation [T7| are easily ob- 
tained for both Rp and Rs. For x^/dof{dof} we obtain 
0.2{4} and 0.4{5} for Rp and Rs respectively. The phys- 
ical results we obtain are: 



l^Pi - IS 
{2S-lS)r 
{3S - lS)r 
{2S - lS)r 



0.820(12) 



= 1.625(39) 



Table XI Most of the errors are obtained directly from 
our fit. The NRQCD systematic error in Rp is taken 



from combining results in Tables VI and VII Since these 
errors are correlated between the numerator and denomi- 
nator of Rp we take the systematic error in i^p to be the 
difference between them. The total NRQCD systematic 
error at each lattice spacing is then included in our fit 
as a correlated error on the data. In fact we find no sig- 
nificant difference whether we include it as a correlated 
or uncorrelated error. We obtain the error in our final 
result from this systematic error by observing the change 
in the final answer from including it or not including it. 
Variation in the NRQCD systematic errors as a function 
of ami) is included in our fit form and the error from this 
estimated from the variation of iii the fit. We use the 
same approach for Rs and take the NRQCD systematic 
error to be the same as for Rp. We might expect some 
further cancellation of errors within Rs because of the 
similarity between the S-wdive states. However, this is 
less true when comparing 35 to 16* than for 2S and IS* 
so we ignore that possibility to be conservative. 

We believe that errors from any mistuning of 7725 are 
completely negligible. Rp and Rs change experimentally 
very little between b and c and we have very well-tuned 
b masses except on the very coarse lattices where our 
mistuning amounts to 4%. 

We also believe that finite volume errors are negligi- 
ble. A study using the heavy quark potential derived 
from a quenched lattice QCD calculation in [36j calcu- 
lated wavefunctions for radially excited T states. None 
of the wavefunctions for the states being considered here 
extended beyond a radius of 1.5 fm and the 2S and IP 
extended little beyond 1.0 fm. When sea quarks are in- 
cluded, as here, the size of the states will be smaller be- 
cause the Coulomb coefficient in the heavy quark poten- 
tial is larger. Thus 1.5 fm is an overestimate for the size 
of the states. The physical extent of our lattices range 
from 2.3 fm for set 1 to 3.8 fm for set 4, so should be large 
enough to contain the T states without any finite- volume 
errors from their being squeezed. 
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In Table XI we include a 0.2% error from electromag- 
netic effects and the possibility of 7^5 annihilation, neither 
of which is included in our calculation. Electromagnetic 
effects we estimated earlier at 1.6 MeV in the IS mass 
and 0.8 MeV (correlated) in the 2S mass. If we take 
the effects on the IP and 3S masses to be much smaller 
then we arrive at a possible error in Rs and Rp of the 
order of 0.1% to 0.2%. 7^5 annihilation affects the spin- 
averaged IS mass, shfting it by approximately 0.5 MeV. 
This amounts to a possible 0.1% effect in Rp, whereas 
Rs is unaffected. 

Our result for the ratio Rp of 0.820(12) is to be com- 
pared with the experimental result 0.8088(23). Agree- 
ment is good within our 1.4% errors. Similarly we obtain 
1.625(39) for Rs to be compared with the experimen- 
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TABLE XI: Complete error budget for the ratios of mass 
splittings, Rp = (l^Pi - 1S)/(2S - 1S)t and Rs = (35' - 
1*S')t/(25' — 1S)t- Errors are given as a percentage of the 
ratio. Errors which are negligible compared to the others are 
indicated by '0'. 



Set amb{ar) 



amsiar) 





Rp 


Rs 


stats/fitting 


1.0 


1.8 


a-dependence 


0.6 


1.2 


m^-dependence 


0.6 


0.5 


NRQCD amb-dependence 


0.1 


0.2 


NRQCD systematics 


0.5 


1.0 


finite volume 








TTib tuning 








electromagnetism/ r]b annihilation 


0.2 


0.2 


Total 


1.4 


2.4 



tal result 1.5896(12). Again agreement is good, but now 
with 2.4% errors, dominated by our statistical/fitting er- 
ror because the 3S state is a doubly excited state. The 
fact that our central value is slightly higher than experi- 
ment for both Rs and Rp is consistent with the expected 
effect of missing terms, included in our errors, as dis- 
cussed in section [Til Dl 

Our result for Rp can be converted to a result for 
M{hb) - M{1S) = 0.461(7) GeV. This can be compared 
to the result 0.440 ± 17+ o^GeV with over double the er- 
ror obtained on configurations including 2+1 flavors of 
sea quarks using the Fermilab heavy quark action [37] . 
The experimental result for M{hb)-M{lS) is 0.4553(17) 
GeV [^Eg. 

Our result for Rs gives M{V') - M(T) = 0.914(23) 
GeV compared to an experimental result of 0.8949(6) 
GeV [35j. We have not included in our error budget any 
effect from coupling of the T^' to virtual decay channels. 
The T^' is 200 MeV below threshold for real decay to a 
pair of B mesons. This is considered large enough for the 
T'^ to be 'gold-plated' and for the decay channel to have 
no impact on the mass. 



2. Tuned b quark masses 

We now return to the tuning of the b quark mass. 
Although not an issue for the mass splittings just dis- 
cussed, it is an important source of systematic error for 
spin-dependent mass splittings. We use our determina- 
tion of the lattice spacing from the T {2S — IS) split- 
ting, given in T able \X\ to convert the kinetic mass values 
given in section |III C| to physical units. As described in 
. Ill C| the appropriate experimental value for com- 



section 

parison is 9.445(2) GeV. 

When this is done we see that the masses are very 
well-tuned except on the very coarse lattices where they 
are 4% high. For small changes in the b quark mass the 
change in kinetic mass is approximately twice the change 
in quark mass, as can be seen from Table [V} For the slight 
changes that we need to make this is a sufficiently good 



1 3.297(11)(35)(7)(16) 

2 3.263(7) (35) (4) (16) 

3 2.696(4)(22)(7)(13) 

4 2.623(7) (22) (7) (13) 

5 1.893(6)(12)(5)(9) 



0.0641(4)(12)(2) 

0.0636(3)(12)(2) 

0.0528(2)(8)(2) 

0.0512(3)(8)(2) 

0.0364(2)(4)(1) 



0.0705(9) (4) 
0.0692(5)(4) 
0.0541(6)(3) 
0.0531(4)(3) 
0.0376(5)(2) 



TABLE XII: Tuned b and s quark masses in lattice units 
on each set of configurations. The second column gives arrih 
and the third arus using the T 2S —IS splitting to determine 
the lattice spacing. The first two errors in these two columns 
come from statistical errors and systematic errors respectively 
in the lattice spacing determination. The third and fourth 
errors in the arrih case are the statistical and systematic errors 
in the determining the kinetic mass. Statistical errors in the 
determination of arris mass are negligible. The third error in 
the arUs case is a correlated 0.3% error from the square of the 
physical value of the rjs mass. The fourth column gives arris 
using the rjs decay constant to determine the lattice spacing. 
The first error is from statistics/fitting and the second is a 
correlated 0.6% error from the square of the physical value of 
the rjs decay constant. 



approximation. We simply adjust the quark mass by one 
half the error in the kinetic mass to obtain the tuned 



quark mass values in lattice units given in Table XII 



Three errors are given in Table XII The first is from 
the statistical error in the kinetic mass determined on 
each ensemble. The second error comes from the total 
error in the determination of the lattice spacing in Ta- 
ble [Xj This includes both statistical and systematic er- 
rors in determining the 26' — 15* splitting. The third error 
is a 0.5% systematic error from NRQCD in the kinetic 
mass obtained from analysis of the dispersion relation in 
section IIII CI 



3. The hyperfine splitting 

The mass difference between the ^Si and ^Sq states is 
an important test of our calculations because it is sta- 
tistically very precise for the ground-state mesons. Con- 
trolling systematic errors is the key issue, and the main 
one is that of radiative corrections to C4, the coefficient 
of the (J • B term in the NRQCD action. At leading order 
the hyperfine splitting is proportional to C4. Our pre- 
vious calculation [2 , with C4 = 1, gave a prediction for 
M(T) - M{r]b) of 61(14) MeV with the error dominated 
by the then-unknown radiative corrections to C4. 

For this calculation we have results for C4 including 
0{as) corrections as well as C4 tuned nonperturbatively. 
These determinations of C4 are described in Appendix [B] 
and Appendix[C]respectively. The values obtained by the 
two methods for C4 are given in Table [XTTTj The nonper- 
turbative values are slightly larger than the perturbative 
ones, but the differences are well within the expectations 
from additional a'^ corrections to the perturbative values 
and/or systematic errors in the nonperturbative values. 
Both sets of values get closer to 1 on the finer lattices, as 
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FIG. 10: Results for the hyperfine splitting, M(T) - M{r]b) 
plotted against the square of the lattice spacing. We show 
results for C4 = 1 (cyan squares) as well as results for C4 set 
equal to its perturbatively improved (red crosses) and nonper- 
turbatively improved values (blue stars). The C4 = 1 results 
include statistical errors only and are shown purely for com- 
parison purposes - they are not included in the fit. The results 
for perturbative and nonperturbative C4 include a correction 
for missing 4-quark operators and mb-mistuning. The errors 
on these points are from statistics, the lattice spacing and the 
tuning of mt. The results for nonperturbative C4 also include 
statistical errors in the determination of C4. Our final physi- 
cal result including our full error budget is given by the grey 
shaded band. The full error budget includes errors from sys- 
tematic uncertainties in setting C4 and from missing terms 
in the NRQCD action. 



TABLE XIII: Values for the coefficient of the a • B term, C4, 
for different lattice spacing values. The error on the pertur- 
bative values is 1 X ttg. The errors on the nonperturbative 
values are statistics, experiment and NRQCD systematics re- 
spectively. We did not extract a nonperturbative value on the 
very coarse lattices. 



Sets 


pert 


^nonpert 


fine 


1.16(5) 


1.18(2)(1)(5) 


coarse 


1.20(7) 


1.28(7)(1)(5) 


very coarse 


1.22(8) 





expected, because they are functions of the strong cou- 
pling constant at a scale related to the inverse of the 
lattice spacing. 



Table XIV gives results for the energies of the T and 
775 for various combinations of values of coefficients in 
the NRQCD action. We also give the mass difference 
between the T and 775 which is the hyperfine splitting. 
This can be more precise than either mass separately be- 



cause we fit both meson correlators together and extract 
the difference directly from the fit taking into account 
the correlations. Where we have fits to a 5 x 5 matrix 
of correlators, as in Tables VIII and |IX[ we give those 



results. In other cases we calculated only a single local 
correlator for each of the T and 7^5 which is quite suf- 
ficient to extract a splitting between the ground state 
masses. Results are not as precise for splittings between 
radially excited states in those cases, however, and we do 
not give them. 

We see from Table |XIV| that changing C4 does have a 
large effect on the hyperfine splitting, approximately in 
line with the expectation of variation as C4. We also see 
that, on sets 3 and 5 where we have data for comparison, 
changing ci^s^e to their 0{as) improved values does in- 
crease the hyperfine splitting slightly. It is a small effect, 
however, of order 2%. Changing the coefficient of the 
Darwin term, C2 also has a small effect of order 1(1)%. 

We conclude from this that we have controlled all of 
the coefficients of terms in our NRQCD action at a 
level required to give few percent errors in the hyper- 
fine splitting from these sources. We have not, however, 
included 4-quark operators in our NRQCD Hamiltonian 
and they can have an impact on the hyperfine splitting 
at an order equivalent to that of ag corrections to C4 [28j . 

In Appendix [C] we give coefficients for the 4-quark op- 
erators and a formula in equation |B14| for the correction 
that they would induce in the hyperfine splitting. Ta- 
ble |XV| gives values for this correction, to be added to 
the results from Table XIV[ on the very coarse, coarse 
and fine ensembles based on using a spin- averaged value 
of the 'wavefunction-at-the-origin', V^(0), for the T and 775 
from our fits. This varies only very little with C4, falling 
by at most 2% from C4 = 1 to our nonperturbative C4 
values, so we ignore this variation. Our correlators are 
normalised by dividing by 6 after summing over 2 spins 
and 3 colors. Then ?/^(0) is given by the amplitude of our 
correlator fits; a(/, 1) for the IS and a(/,2) for the 2S 
from equation [71 

We see from Table [XVl that the correction is substan- 
tial on the very coarse lattices and very small on the fine 
lattices, because of the variation of the coefficients di and 
d2 with ami). The corrected results are shown in Fig- 

1. We 



ure 



10 along with the uncorrected results for C4 
see that there is a substantial difference between the re- 
sults coming from the change in C4 and, to a lesser extent, 
from the correction for the 4-quark operator. The strong 
dependence on the lattice spacing seen in the C4 = 1 
results is reduced, in line with the expectation that im- 
proving an effective theory should reduce the cutoff de- 
pendence. 



An additional small factor in Figure 10 is that we have 
corrected results for slight mistuning of the b quark mass 
and we have included the error from the quark mass tun- 
ing in the hyperfine splitting error. The hyperfine split- 
ting is expected to be approximately inversely propor- 



tional to the quark mass and this is seen in Table XIV 



We assume this relationship to make small adjustments 
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Set 


airib 


Cl,5,6 


C2 C4 




aEr 


aEr — ciEr^^ 


1 


3.42 


as 


1 1 


0.25080(5) 


0.28532(6) 


0.03452(8) 


1 


3.42 


as 


1 1.22 


0.21432(5) 


0.26400(6) 


0.04968(7) 


1 


3.5 


as 


1 1 


0.25015(6) 


0.28392(9) 


0.03377(10) 


2 


3.39 


as 


1 1 


0.25361(3) 


0.28809(3) 


0.03448(4) 


2 


3.42 


as 


1 1 


0.25344(5) 


0.28759(5) 


0.03416(6) 


3 


2.66 


1 


1 1 


0.25529(4) 


0.28626(6) 


0.03097(7) 


3 


2.66 


as 


1 1 


0.26096(3) 


0.29245(3) 


0.03149(4) 


3 


2.66 


as 


1.25 1 


0.25627(24) 0.28728(33) 0.03101(41) 


3 


2.66 


as 


1 1.25 


0.20943(3) 


0.25628(3) 


0.04684(5) 


3 


2.66 


as 


1 1.20 


0.22040(5) 


0.26394(7) 


0.04354(4) 


3 


2.68 


as 


1 1 


0.26108(7) 


0.29249(9) 


0.03141(11) 


3 


2.7 


1 


1 1 


0.24375(8) 


0.27483(11) 


0.03108(13) 


4 


2.62 


as 


1 1 


0.26524(2) 


0.29681(2) 


0.03157(3) 


4 


2.62 


as 


1 1.25 


0.21289(2) 


0.25978(2) 


0.04689(2) 


4 


2.66 


as 


1 1 


0.26546(3) 


0.29662(4) 


0.03116(5) 


5 


1.91 


1 


1 1 


0.24652(3) 


0.27153(5) 


0.02501(6) 


5 


1.91 


as 


1 1 


0.25851(2) 


0.28405(2) 


0.02554(3) 


5 


1.91 


as 


1 1.10 


0.23204(2) 


0.26206(3) 


0.03003(4) 


5 


1.91 


as 


1 1.15* 0.21772(22) 0.24984(40) 0.03213(30) 


5 


1.91 


as 


1 1.16 


0.21519(2) 


0.24802(4) 


0.03283(2) 


5 


2.0 


as 


1 1 


0.25935(3) 


0.28397(4) 


0.02462(5) 



TABLE XIV: Fitted energies for ground state T and 775 mesons on all configuration sets. The column ci,5,6 denotes whether 
the 0{as) improved coefficients were used in the action. Various values of C2 and C4 have also been used as indicated. C3 = 1 
except for the case indicated by * in which C3 = 0.96. Where possible the result from the full 5x5 matrix fit was taken. 
Otherwise, the values from the kinetic mass fits were used. In those cases a smaller number of configurations and/or time 
sources was sometimes used and this is reflected in the statistical errors. 



TABLE XV: Corrections to the IS and 25' hyperfine split- 
tings from spin-dependent 4-quark operators missing from our 
NRQCD action. We use equation |B14| inserting values for 
'^(O) from our fitted results and values for av(7T/a) from Ta- 
ble |XXII[ We convert to physical units using lattice spacing 
values from Table [Xl 



ing to the second and third errors in Table XIII and re- 



Sets 


Correction to 15 


Correction to 25 




hyperfine (MeV) 


hyperfine (MeV) 


fine 


-1.7 


-1.0 


coarse 


5.2 


3.4 


very coarse 


12.9 


8.3 



based on the tuned b masses from Table XII The largest 
effect is a 4% one on the very coarse lattices. The lattice 
spacing error from the quark mass tuning is correlated 
with the lattice spacing error on the hyperfine splitting 
because of this inverse relationship. The lattice spacing 
error therefore appears with a factor of 2 in the hyperfine 
splitting. 

To obtain a physical result for the hyperfine splitting 
we then combine results with perturbative and nonper- 
turbative values of C4 allowing for systematic differences 
between them from uncertainties in the determination 
of C4. We must also allow for uncertainties from higher- 
order 4-quark operator effects and for lattice spacing and 
sea quark mass dependence. The nonperturbative C4 re- 
sults are given a correlated systematic error correspond- 



membering that the hyperfine splitting is related to C4. 
Similarly the results for perturbative C4 are given a sepa- 
rate correlated systematic error corresponding to the 



errors given in Table XIII We allow for higher order 4- 
quark operator effects with a correlated systematic error 
of size 6af |?/^(0)p/m^ with a coefficient of possible size 
±1 ± ln(am5). This does not assume that the small co- 
efficient seen at 0{al) on the fine lattices is repeated at 
higher order. 

We allow for lattice spacing and sea quark mass effects 
as in the fit function of equation 17 The prior on sea 
quark mass effects is now taken to allow 15% effects for 
mi ~ rus. This reflects the fact that a 40% difference 
was seen between quenched and dynamical results in [2] . 
In fact sea quark mass effects in our data are small. Fig- 
ure [TT] shows a comparison of the hyperfine splitting as a 
function of the light sea quark mass for the case C4 = 1 
where we have a complete set of data. Although these 
results are not used to determine our final answer for 
the hyperfine splitting they do provide a useful compar- 
ison between ensembles at the same lattice spacing and 
different light quark mass since the effects of C4 are in- 
dependent of sea quark mass. The results are adjusted 
for b quark mass mistuning and include errors from the b 
quark mass and the lattice spacing. Variation with light 
quark mass is at most 2 MeV. Note that we are using 
much lighter sea quark masses than in previous calcu- 
lations [2 [T^; indeed mi is within a factor of 3 of its 
physical value. 
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FIG. 11: Results for the hyperfine splitting obtained for C4 = 
1 (not used in our fit for the physical hyperfine splitting) 
compared as a function of sea light quark mass in units of 
the strange quark mass. Results are given for two values of 
mi/ms on very coarse (sets 1 and 2) and coarse (sets 3 and 
4) lattices. The errors on the points include statistical/fitting 
errors, lattice spacing errors and mt tuning errors. 



TABLE XVI: Complete error budget for the 15* hyperfine 
splitting and the ratio of the 2S to the 15* hyperfine splittings. 
Errors are given as a percentage of the final result. * Note 
that the mt tuning uncertainty does not include the lattice 
spacing uncertainty in m^. Since that is correlated with the a 
uncertainty on converting the hyperfine splitting from lattice 
to physical units, they must be handled together and both 
are included in the a- uncertainty. 





M(T) - M{r]b) 


Rh 


stats/fitting 


0.1 


4 


a-dependence 


1.5 


5 


a-uncertainty 


0.5 





771; -dependence 


3 


3.5 


NRQCD amb-dependence 


2 


0.5 


NRQCD 


10 


5 


NRQCD C4 uncertainty 


7 





NRQCD 4-quark uncertainty 


2 


1 


rub tuning* 


0.1 





r]b annihilation 


1 


0.5 


Total 


13 


9 



We obtain a physical value for the hyperfine splitting 
from the above fit of 70(6) MeV. Fitting the results for 
perturbative C4 on their own gives a consistent 67(7) MeV 
and the results for nonperturbative C4 alone gives 75(9) 
MeV. Our best result therefore comes from combining 
the two. An additional 10% error must be allowed for 
higher order (v^) spin-dependent terms in the NRQCD 
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2011 



HPQCD (NRQCD) 
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Meinel (NRQCD) 
2010 [12] 



Fermilab/MILC 
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Hyperfine splitting (MeV) 

FIG. 12: Comparison of results for the hyperfine splittings, 
M(T) - M{r]b) and M{V) - M{r]'b), from different full lattice 
QCD calculations. Filled symbols indicate the IS hyperfine 
and open symbols, the 2S hyperfine. Circles indicate pre- 
dictions and triangles postdictions. The top (red) points are 
the new results from this paper, and the points below that 
(pink) are from 2 , when the IS hyperfine was a prediction. 
The third line gives results from |12j . Two results are given 
for the 2S hyperfine; that from a ratio to the IS hyperfine, 
as here, and that from a ratio to the combination of P-wave 
spin splittings sensitive to C4 (see Appendix[C]). The top three 
results use the NRQCD formalism for b quarks; the bottom 
(cyan) result uses the Fermilab heavy quark action 33 • The 
black dashed lines mark the current experimental average [35] . 



action, giving a final result of: 

M(T) - M{r]b) = 70(9)MeV. 
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Our complete error budget is given in Table XVI [ The 
shift for the effect of 775 annihilation is included in our 4- 
quark operator correction (as discussed in Appendix [C|) 
but we separate out an error for that from the rest of the 
4-quark operator error. 



In Figure 12 we compare our new result for the hy- 
perfine splitting to earlier full lattice QCD results and 
to experiment [35]. Earlier results using NRQCD are: 
61(14) MeV from [2l with a treelevel NRQCD ac- 
tion and 60.3(7.7) MeV from [12^ using an NRQCD 
with spin-dependent terms and C4 determined from P- 
wave splittings but with no 4-quark operator corrections, 
which are potentially more significant than terms, or 
errors from them. The result obtained from the Fermilab 
heavy quark action [37 is 54.0(12.4) MeV. For this ac- 
tion the hyperfine splitting is sensitive to the coefficient 
of the 0{a) improvement term known as the clover term. 
In principle this coefficient does not have to be tuned 
but the approach to the continuum limit is slow. Here 
it was taken to have its tree-level value after tadpole- 
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improvement and the quark mass was tuned using the 
spin- averaged kinetic mass of the Bs and B* mesons. 

Ah four lattice QCD results agree well with the current 
experimental average of 69.3(2.8) MeV [35^ obtained from 
averaging results from experiments on radiative transi- 
tions to ?75 from and T^' [39H4T] . Preliminary exper- 
imental results using radiative transitions from the 
indicate a somewhat lower value [42 . 

Our new result above contains the most complete anal- 
ysis at 0{v^) in the NRQCD action. To improve it 
would require the inclusion of spin-dependent operators 
at 0{v^). The effect of spin-dependent operators was 
studied in [12] , taking ratios of the hyperfine splitting to 
P-wave spin splittings to cancel the effect of C4. Those 
results indicate that spin-dependent terms tend to re- 
duce the hyperfine splitting by about 10%. We have in- 
cluded a (symmetric) 10% error in our results to account 
for missing terms. 

The hyperfine splitting has also been calculated us- 
ing continuum QCD perturbation theory [43 . A con- 
siderably smaller result is obtained of 41(14) MeV. This 
is not in disagreement with the nonperturbative lattice 
QCD results given the size of the errors. It has been 
suggested, however, that the inclusion of radiative cor- 
rections to C4 in the lattice NRQCD calculation would 
reduce the value of the hyperfine splitting obtained [44j. 
That expectation was based on an incorrect analysis of 
the form of C4 in the lattice NRQCD calculation and we 
see indeed from the results given here that the inclusion 
of radiative corrections to C4 has had the opposite effect 
and increased our value for the hyperfine splitting. 

The best way to study the 25* hyperfine splitting, 
M(T') - M{r]l) is through the ratio to the 16' hyperfine 
splitting. We define: 



Rh 



M{V)-M{i,) 
M(T) - M{r],) ' 
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Then C4 effects cancel as can be seen in the numbers in 
Tables |VIII| and |IX] and plotted in Figure [13] In the 
Figure we have corrected the results for missing 4-quark 
operator effects which are slightly different in the IS and 
2S states and so have an effect on the ratio. This is 
at most 4%, on the very coarse lattices. We have made 
no correction for the slight mistunings of m5 since they 
should largely cancel in this ratio. 

Again we extract a physical result for Rh by allowing 
for mi and a dependence as in equation [17] Here we 
relaxed the priors on the a-dependence so that they all 
had the form 0.0 ± 1.0. The prior on m^-dependence 
(in units of m^) was taken as 0.0 ± 0.15 as for the IS 
hyperfine, allowing a 15% change from mi = ms down to 
the physical point. 

Our final physical value is: 



Rh = 0.499(42). 
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The full error budget is given in Table |XVI| where we 
allow a 5% error for missing spin-dependent terms. 
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FIG. 13: The ratio of hyperfine splittings for 2S and 15* 
states, (M(T') - M(ryb))/(M(T) - M{r]b)), plotted against 
the square of the lattice spacing. Points have been corrected 
for missing 4-quark operator effects. The shaded band shows 
our final physical result including our full error budget. 



allowing for some cancellation of effects between the 
IS and 25* hyperfine splittings. 

Combining our result for Rh with the current experi- 
mental average for the IS hyperfine splitting gives a re- 
sult for the 2S hyperfine splitting of 35(3) (1) MeV, where 
the second error comes from the experimental IS split- 
ting. We predict the mass for the 7^[, to be 9988(3) GeV. 



Figure [12] compares our result for the 2S hyperfine 
splitting to earlier predictions. Meinel [12] gives 23.5(4.7) 
MeV from a ratio to P-wave spin splittings and 28.0(4.7) 
from a ratio to the 15* hyperfine splitting, as used here. 
He includes the effect of spin-dependent operators but 
without such a complete analysis as we have done here 
of operators. The conclusion from Figure 12 is that 



operators. The conclusion from Figure 
lattice QCD results give a fairly clear prediction of the 
2S hyperfine splitting around 30 MeV, half the result for 
the IS hyperfine splitting. 

Figure [13] compares our result for the bottomonium ra- 
tio of 2S to IS hyperfine splittings to the experimental 
charmonium ratio of 0.421(35). Our bottomonium result 
is somewhat higher, but not in disagreement with this 
value. This indicates that the heavy quark mass depen- 
dence in the IS and 2S hyperfine splittings is very similar 
over the wide range of quark masses from c to 6 so that 
the ratio remains the same. 

The 6'-wave hyperfine splittings can be compared to 
the much smaller result for the P-wave states. In Ap- 
pendix [C] we determine the P-wave hyperfine splitting to 
be 2(2) MeV, consistent with zero. 
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IV. THE r]s MASS AND DECAY CONSTANT 

To complement the computation in the previous sec- 
tion, the lattice spacing was also determined using the 
decay constant of the fictitious rjg meson, /^^. This is a 
pseudoscalar particle consisting of an ss pair whose prop- 
erties can easily be computed in lattice QCD. It is par- 
ticularly suitable for fixing the lattice spacing since there 
are no u/d valence quarks meaning that the error coming 
from the chiral extrapolation to physical u/d masses is 
small. The "physical" values of the M^^ and /^^ have 
to be fixed by comparison to M^^, M^, /^^ and fx as 
in [21 . This requires a simultaneous chiral and contin- 
uum extrapolation for the masses and decay constants of 
the TT, K and 7^5. Previously we found, on ensembles in- 
cluding 2+1 fiavors of sea quarks, that properties of the 
r]s were very close to those expected from leading order 
chiral perturbation theory, i.e. M^^ ^ 2M|^ - and 
/^^ ^ 2fK — fn- We re-examine that issue here on these 
ensembles containing 2+1+1 sea quarks. 



A. Simulation details and Fitting 

s and u/d valence quark propagators were calculated 
on the ensembles given in Table [l] using the same HISQ 
action as used for the sea quarks. The HISQ action [16 
is a further improved version of the improved staggered 
(asqtad) action that reduces discretisation errors coming 
from staggered taste effects for that action by about a 
factor of 3. The improved staggered action smears the 
gluon fields that appear in the quark action in a very spe- 
cific way to reduce high-momentum taste-exchange inter- 
actions but without increasing discretisation errors [10]. 
The HISQ action takes this one step further by perform- 
ing two smearing steps. The original version of the HISQ 
action used an SU(3) projection of the smeared links be- 
tween the two smearing steps. However, this caused dif- 
ficulties for the updating algorithm when these quarks 
were included as sea quarks [19] . Instead the sea quarks 
here use the HISQ action but with only a U(3) projec- 
tion between the smearing steps. Whether U(3) or SU(3) 
projection it makes very little difference to the spectrum 
of mesons from HISQ quarks and so all the good features 
of the HISQ action demonstrated in [16 remain essen- 
tially unaltered. Note that the HISQ action does not 
use tadpole-improvement - the U(3) (or SU(3)) projec- 
tion effectively takes care of large tadpole contributions 
in the same way that the use of the uq parameter does 
in the NRQCD and gluon actions. 

The parameters of the valence quark propagators are 



listed in Table |XVII| We took the light quark mass to 
be the same as that in the sea (except for a small differ- 
ence on set 3), but we retuned the valence strange quark 
masses slightly to allow for mistuning of the strange sea 
quark mass (of course the final well-tuned s quark mass 
values cannot be decided until a value for the lattice spac- 
ing is determined so we will revisit this issue at the end of 



TABLE XVU: Valence light and strange quark mass parame- 
ters on each ensemble. The valence light quark masses are the 
same as in the sea (given in Table except for a slight dif- 
ference on set 3. The valence strange quark masses have been 
retuned slightly to be closer to the physical values. Columns 
4 and 5 give the number of configurations used from each en- 
semble and the number of time sources for propagators per 
configuration. 



Set 


amr^ 






nt 


1 


0.013 


0.0688 


1021 


16 


2 


0.0064 


0.0679 


1000 


16 


3 


0.01044 


0.0522 


1053 


16 


4 


0.00507 


0.0505 


1000 


16 


5 


0.0074 


0.0364 


1008 


16 



this section). We used delta function random wall sources 
as for the 1-smeared b quark propagators discussed in sec- 
tion 



IHl We also used 16 evenly-spaced time sources per 
configuration to increase statistics. The starting posi- 
tion of these time sources was shifted from configuration 
to configuration in the ensemble. 

The light meson pseudoscalar correlators are calcu- 
lated by combining the light and strange quark propaga- 
tors. Here we use the goldstone mesons, made with the 
local 75 operator. Then the correlators are simply given 
by the squared modulus of the propagators summed over 
a time-slice to project onto zero momentum. The correla- 
tors were binned over all time sources on a configuration. 
To study the autocorrelations between configurations 



we proceed as in subsection |III B|t o calculate the autocor- 
relation function Cat- Figure] 14] shows Cat against AT 
for both the tt and 77^ correlators at a source-sink lattice 
time separation appropriate to our fits. This time separa- 
tion is increased as the lattice spacing decreases to remain 
approximately physical. We see from the Figure that the 
r]s shows little autocorrelation, although more than was 
visible for the T in Figure [ij We expect longer autocorre- 
lation times for lighter mesons because they have a larger 
spatial extent and therefore decorrelation of the gluon 
field configuration on relevant spatial scales takes longer 
in Monte Carlo time. The tt correlators on the coarse 
and very coarse lattices show a similar autocorrelation 
function to the rjs. However, on the fine lattices where 
autocorrelations might be expected to be worst, there are 
clear signs that neighbouring configurations in time are 
correlated. Note the difference between our result and 
that of [19j. There very little autocorrelation was seen in 
the TT meson correlator, but fewer time sources were used 
per configuration (typically 4). Both here and in [19 the 
time sources are moved randomly from one configuration 
to the next. However, with 16 time sources there is not 
much scope for a large shift in time between configura- 
tions. To reduce autocorrelations we bin all of our fine 
light meson correlators by a factor of 8 in configuration 



number before fitting. From Figure 14 this can be seen 
to reduce the autocorrelation function well below e~^. 
The fitting method used for the correlators was the 



23 




-0.2 ' ' ' ' ' ' -0.2 ' ' ' ' ' ' 

02468 10 02468 10 

AT AT 

FIG. 14: Autocorrelation function Cat for tt (left) and r/s (right) correlators. These are made with a S function random wall as 
described in the text. The key to results from different ensembles is (color online): set 1, green plus; set 2, orange cross; set 3 
blue star; set 4 pink open square; set 5, red open circle. The correlators are evaluated at lattice time separation t/a = 6 on very 
coarse lattices (sets 1 and 2), t/a = 8 on coarse lattices (sets 3 and 4) and t/a = 10 on fine lattices (set 5). This corresponds 
to a t value, approximately constant in physical units across the lattice spacing values, where the tt and r]s correlators have 
reached the ground-state plateau. AT is given in units of configuration number in the ordered list for each ensemble. 



same as for the Upsilon correlators but with only a single 
source and sink smearing. For meson correlators made 
from relativist ic quarks the fit function takes a 'cosh' 
form rather than simple exponentials because of prop- 
agation in both time directions. For staggered quarks 
in general we have to include an additional oscillating 
term from opposite parity mesons that couple through 
the time-doubler quark. The fit function then becomes: 



from the fit using 



G„{t) = ^afc(e-^^* + e-^''(^-*)) (21) 

fc=0 



ko=0 



The oscillating piece is absent for the tt and rjs because 
the valence quark and antiquark have equal mass and 
the oscillation cancels. It is necessary to include it for 
the K meson. For each ensemble a simultaneous fit to 
all three correlators was performed using the appropriate 
form for each. This allowed us to take into account the 
correlations between the fit results for each meson in our 
subsequent chiral extrapolations. We use the full range 
of t values in the fit apart from the first 3-5 time-slices. 
Priors for energies and amplitudes are chosen as for the 
T fits described in section IIIIBI We take results from 
4 exponential fits, where ground-state masses and their 
errors have clearly stabilised. 

The results that we use are the ground-state meson 
masses and amplitudes, i.e. /c = in equation [21] The 



meson masses are given by the parameter Eq from each 
fit, since there is no energy offset for staggered quarks as 
there is for NRQCD. The decay constants are extracted 




(22) 



for a meson containing quarks a and b with masses 
rria^mi)^ ground state mass Eq and ground state ampli- 
tude, ao, from the fit form above. Equation [22] uses the 
PCAC relation, valid for staggered quarks, to relate the 
matrix element of the pseudoscalar density to that of 
the temporal axial current and therefore the decay con- 
stant. The existence of the PCAC relation means that 
the temporal axial current is absolutely normalised and 
there is no uncertainty from lattice to continuum current 
matching factors as there can be in some other quark 
formalisms. 



B. Results and chiral extrapolations 



Our results for the tt, and r]s meson masses and de- 



cay constants are listed in tables XVIII and XIX We also 
give various ratios that are useful indicators of the sensi- 
tivity of the rjs parameters to the chiral extrapolation in 
the u/d sea quark mass, and to the lattice spacing. 

We fit the three decay constants and meson masses 
simultaneously using SU(3) chiral perturbation theory, 
adapted to include discretisation effects. The usual ap- 
proach is to use values for the decay constant and me- 
son mass in GeV, having chosen a value of the lattice 
spacing on each ensemble. Extrapolation to the point 
where and Mk take their physical values then allows 
comparison to experiment of the resulting values for 
and fx- Here instead we use values for ri/a to fix the 
relative lattice spacing between ensembles and keep the 
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oet 


1 


2 3 


4 


5 


aM^ 
aMx 


0.23637(15) 
0.41195(17) 
0.53361(14) 


0.16615(7) 0.19153(9) 
0.39082(9) 0.32781(10) 
0.52797(8) 0.42351(9) 


0.13413(5) 
0.30757(7) 
0.41476(6) 


0.14070(9) 
0.23933(11) 
0.30884(11) 


Ml/{2M-t, 


-Mi) 1.00426(43) 1.00317(28) 1.00636(34) 1.00474(26) 1.00660(27) 



TABLE XVIII: Values for the ground state masses in lattice units (Eq from eq. 21) for tt, K and r/s mesons. The fourth 
row gives the ratio of the square of the r]s mass to a combination of tt and K masses that would be 1 in leading order chiral 
perturbation theory. 



Set 1 2 3 4 5 

aU 0.11183(9) 0.10511(5) 0.09075(5) 0.08451(4) 0.06621(5) 

afK 0.12689(8) 0.12268(4) 0.10185(5) 0.09788(3) 0.07427(4) 

afru 0.14199(6) 0.14026(3) 0.11312(4) 0.11119(2) 0.08238(4) 

fK/fn 1.13467(58) 1.16717(38) 1.12231(38) 1.15819(35) 1.12170(39) 

fr^JU 1.26974(80) 1.33442(59) 1.24653(68) 1.31568(53) 1.24416(69) 

fvs/i'^fK - f^) 1.00031(62) 1.00007(27) 1.00154(45) 0.99948(26) 1.00061(32) 

fr^jMr^^ 0.26609(11) 0.26566(6) 0.26711(10) 0.26809(6) 0.26674(12) 



TABLE XIX: Values for the ground state decay constants in lattice units (derived from ao in eq. 21 as described in the text) 
for TT, K and rfs mesons. We also give various ratios of decay constants obtained from the simultaneous fit. The sixth row 
gives the ratio of the rjs decay constant to a combination of tt and K decay constants that would be 1 in leading order chiral 
perturbation theory. 
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FIG. 15: The pseudoscalar decay constants plotted against 
the ratio of squared pseudoscalar masses that is approxi- 
mately equal to mi/rris. The points have been adjusted for 
finite volume effects and for mistuning of the strange quark 
mass. The lines are from the tuned fit function at each lat- 
tice spacing, with results increasing in value from very coarse 
(blue) to fine (red) (color online). The top (black) line is 
the a = curve and the black leftmost data points give the 
experimental value for f^^ and fx given values for Vud and 
Vus [35]. 



physical value of ri as a parameter to be obtained from 
the fit. The value for ri is determined by the require- 
ment to match /^^ and fx from experiment in the chiral 
limit where the experimental values are included as extra 
pieces of 'data' for the fit. The experimental values for 
/t^ and Jk come from experimental measurement of the 
leptonic decay rate and values of Vud and Vus taken from 




0.1 0.2 0.3 0.4 0.5 

Ml/{2Ml-Ml) 

FIG. 16: The ratio of fn^ to 2/k — /tt, which would be 1 in 
leading order chiral perturbation theory. The ratio of squared 
meson masses on the x-axis corresponds approximately to 
mi/ Tils- The blue, green and red points and fit curves cor- 
respond to very coarse, coarse and fine lattices respectively 
(color online). The black line is the continuum, a = 0, fit 



elsewhere. We use 



= 0.1304(2)GeV 
fK = 0.1561(9)GeV. 



(23) 



The meson mass values that go with these decay con- 
stants in a world appropriate to lattice QCD without 
electromagnetism and in which ruu = md are [45^ : 

I (M^o + M|,+ - (1 + Ab)(M2+ - M^o)) . 



(24) 



K 



25 



^ 1.02 
t 1-00 




0.1 0.2 0.3 0.4 0.5 

Ml/{2Ml-Ml) 



FIG. 17: The ratio of to 2Mi - M;, which would 

be 1 in leading order chiral perturbation theory. The ratio 
of squared meson masses on the x-axis corresponds approx- 
imately to mi/rris. The blue, green and red points and fit 
curves correspond to very coarse, coarse and fine lattices re- 
spectively (color online). The black line is the continuum, 
a = 0, fit curve. 



We take A^;, which parameterizes the violations of 
Dashen's theorem, to have the value 1 ± 1. The decay 
constants are already defined to be results in pure QCD, 
provided electromagnetic effects have been removed from 
the experimental leptonic decay rates [35 . The residual 
error in Jk from the fact that it is the decay constant of 
the whereas the K mass in equation 25 is the isospin 
average is less than 0.1% [45 , so we ignore it here. 



Our fit also returns a physical value for /^^ and M^^ . 
This can be used in subsequent analyses to tune the s 
quark mass and to fix the lattice spacing. It is a good 
meson to use for this purpose because its parameters are 
very insensitive to the sea quark masses, as we see in 
Tables IXVIIII and IXIXl 

The analysis is the same as that used in [21 except for 
two improvements. The first is that the fit is simplified 
because HISQ quarks are used in both the valence and 
sea sectors. The second is that we include correlations 
between all of the decay constants and meson masses on 
each ensemble by feeding into the fit the covariance ma- 
trix that resulted from the simultaneous fit to all three 
meson correlation functions. Below we provide a brief 
description of the chiral/continuum extrapolations fol- 
lowing [21]. 

On each ensemble the decay constants and meson 
masses are a function of the masses of the valence quarks 
for that meson and of the masses of the sea quarks for 
that ensemble, with the coefficients of the mass depen- 
dence constrained to be the same for the tt, K and rjs 
mesons. Because the quark masses run with energy scale 
it simplifies the fits to use a dependent variable related 
to the square of appropriate goldstone meson masses in- 



stead the quark mass. Thus we write 

Ml/2 



xi 



A2 



(25) 



where provides the cut-off scale of the chiral expan- 
sion, 



A, 



47r/,/V2. 



Similarly 



Ml 



Mi/2 



(26) 



(27) 



In the cases where the sea and valence quark masses dif- 
fer, the X parameters for the sea quark masses are ob- 
tained from those of the valence masses by rescaling in 
proportion to the quark mass. Then the decay constant 
made from valence quarks a and h takes the functional 
form 



f(rr ^sea sea 



f 



NLO 



/lat7 



(28) 



where f^^^ is the full partially quenched chiral perturba- 
tion theory formula at next-to-leading order [46] and Sf^ 
and ^/lat include possible correction terms coming from 
higher order terms in the quark mass and finite lattice 
spacing corrections. Each of the terms contains a set of 
unknown coefficients which are given prior constraints in 
our fit allowing us to test their effect on our final result. 

yNLO includes terms proportional to the squares 
of the appropriate meson masses as well as logarithmic 
terms that appear in combinations such as, for exam- 



ple, {Xa 



^). The logarithmic terms 



are corrected for finite volume effects through the use of 
finite volume chiral perturbation theory. The finite vol- 
ume correction is significant for /t^, particularly on set 
1 where M^L = 3.8 and the finite volume correction is 
1.8%. For the other sets, with M^L > 4, the correction 
ranges from 0.4% to 0.7%. For fx and fr^^ the correction 
is much smaller. It is at its largest on set 1 with 0.7% for 
/xand 0.2% for /^^. 

Sf^ includes polynomial dependence on various combi- 
nations of the Xi up to and including xj terms [21j. Most 
of these terms only matter for the s quark and they allow 
for differences between the s and / sectors within SU(3) 
chiral perturbation theory. Since Xg = 0.17 including 
xj terms means that missing terms at x^ are O(10~^), 
smaller than our statistical errors. It is sufficient to in- 
clude polynomials because we cannot distinguish high- 
order logarithms from polynomials over this range in Xi. 

J/iat allows for dependence on powers of the square of 
the lattice spacing, since this is the form that discretisa- 
tion errors take for staggered quarks. We include terms 
up to (aAqcD)^ where Aqcd is taken to be (9(0.6GeV). 
The coefficients of the a-dependence are also allowed to 
have dependence on valence and sea mass dependence. 
This includes dependence on log(x^) to model discreti- 
sation errors coming from staggered taste-changing ef- 
fects 2L. 
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TABLE XX: Complete error budget for ri, /^^, M^^ and 
fr]sl^r]s- Errors are given as a percentage of the physical 
value. Errors which are negligible compared to the others are 
indicated by '0'. 



/^^ and M^^ are 





r\ 








stats/fitting 


0.24 


0.16 


0.07 


0.18 


a-extrapolation 


0.46 


0.14 


0.03 


0.16 


-extrapolation 


0.09 


0.12 


0.04 


0.11 


finite volume 


0.04 











ri/a 


0.73 


0.12 


0.02 


0.12 


initial r\ uncertainty 


0.26 


0.02 





0.02 







0.05 


0.14 


0.09 


Total 


0.90 


0.28 


0.17 


0.30 



The terms in the chiral expansion are generally written 
so that the coefficients are expected to be (9(1). For these 
coefficients we take the prior in our fit to be ± 1. This 
is true for the higher order terms in the chiral expansion 
that relate to a-dependence and mass-dependence, except 
where the masses involved are sea-quark masses and then 
the prior is taken as ± 0.3, simply because sea-quark 
effects are typically suppressed over valence quark effects 
by this amount. The prior on the bare decay constant 
parameter in chiral perturbation theory, /o, is taken as 
0.11 ± 0.02. In fact the parameter that is tuned by the 
fit is log(/o) in order to keep /o positive. The prior on 
log(/o) is then taken as —2.2 ± 0.18. The priors for the 
coefficients ^4,5,6,8 that multiply analytic terms at NLO 
in chiral perturbation theory are taken as ± 0.01. 

The meson masses are fitted simultaneously with the 
decay constants feeding in the 6x6 covariance matrix on 
each ensemble. The leading behaviour in chiral perturba- 
tion theory for the meson masses is now trivial. However 
the chiral fit, which shares some of the same coefficients 
as that of the decay constants [46,, allows us to fix the 
higher order behaviour as a function of sea and valence 
masses. In particular it allows us to fix the behaviour 
of the r]s mass as the tt and K masses vary, so that we 
can obtain its value at the physical point. The priors in 
the chiral fit for the meson masses take the same form as 
described above for the decay constant. 

The fitting forms above were extensively tested for ro- 
bustness against both real and fake data in [17] and [21] . 

The 



The results of our fit are shown in Figure 15 
data points, adjusted for finite volume effects and for 
the slight mistuning of the valence and sea strange quark 
masses, are plotted as a function of xi/xs. The fit lines 
at each value of the lattice spacing are shown along with 
the a = line. At the physical value for xi/xs we 
give the experimental values for /^^ and /k- This plot 
should be compared with Figure 4 in [21 . It is evident 
that these 'second generation' configurations have signif- 
icantly smaller discretisation errors [19]. 

The fit has a x^/dof value of 0.3 for 36 degrees of 
freedom. The fitted value of /o is exp(-2.174 ± 0.028), 
in agreement with SU(3) chiral fits using asqtad improved 
staggered quarks [47 . The resulting physical values for 



0.1819(5)GeV 
0.6893(12)GeV 



fvs/Mr^s = 0.2638(8). 



(29) 



These are in agreement with the results obtained on 
n J = 2 + 1 dynamical asqtad configurations [21] but con- 
siderably more accurate because our statistical precision 
is improved, and we have smaller continuum and chiral 
extrapolation errors. These last two are reduced because 
of the improvements in the gluon field configurations and 
because we are working closer to the chiral limit. Com- 
plete error budgets for /^^ , M^^ and their ratio are given 



in Table IXXl 

The results in equation |29| are very close, but in fact 
differ significantly from the expected result from leading 
order chiral perturbation theory. This is illustrated in 
Figures 16 and 17 in which the ratios frjJi^fx — /tt) and 
M^J{2M]^ — Mj) are plotted against xi/xg. Both ratios 
are very flat in xi/xs^ never differing by as much as 1% 
from 1. We determine the physical values for the ratios 
to differ significantly from 1, however, with results: 

UJ(2fK-U) = 0.9977(6) 
MlJ(2Mj^- Ml) = 1.0070(18). (30) 

Since our fit uses ri/ato set the relative lattice spacing 
we can determine a value for ri from the final match with 
experiment for and fx- We obtain 



ri(/^J = 0.3209(29)fm. 



(31) 



The error budget for ri is given in Table XX This phys- 
ical result for ri agrees with the value obtained from 
the same analysis on n/ = 2 + 1 dynamical asqtad lat- 
tices [21], but is almost twice as accurate. 

We can use /^^ and M^^ to determine the s quark mass 
and lattice spacing on each ensemble. This is done by 
tuning the s quark mass so that fr^jMr^^ takes the value 
in equation |29| and then the lattice spacing is read from 
the value of fr^^ . We do this here retrospectively by using 
our chiral fits to tune the s quark mass and work out the 
corresponding changes in /^^ and M^^ . For simplicity the 
sea light quark masses were also retuned to the physical 
value. The values of a obtained are given in Table [Xj 



In Table |XII| we give tuned values of arris on each en- 
semble as a result of tuning the rjs to the physical value 



given in equation 29 Over the short range needed for the 
retuning the relationship ex M^^ works very well. We 
give results for both the case of using the rjs decay con- 
stant to fix the lattice spacing and of using the T 2S — IS 
splitting. The values of arris obtained from the two meth- 
ods differ substantially on the very coarse lattices but 
come into agreement on the fine lattices as expected. 

The errors in the tuned values of arris are dominated by 
the errors in the lattice spacing. The relative error in a is 
doubled in arris because the quark mass is proportional 
to the square of the meson mass. When the quark mass 
is converted to physical units one factor of the lattice 
spacing error disappears. 
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To extract a physical value for ri using the T results 
we use the same functional form for the fit as was used 
earlier for Rs and Rp^ equation [T7| This includes an 
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FIG. 18: Values for the heavy quark potential parameter n 
obtained by combining values for ri/a from MILC with either 
of our two methods for determining the lattice spacing. The 
red plus symbols correspond to using the r/s and the blue stars 
to using the T 2S — IS splitting (these points do not include 
the NRQCD systematic error which is correlated between the 
points). The black line with light red error band corresponds 
to the final value for ri from the combined fit to the results 
from both methods and includes the total error. 



V. n 

The values of the heavy quark potential parameter, ri, 
can be determined by combining the values for ri/a from 
MILC given in Table [T] with the values for the lattice spac- 
ing given from our two different methods in Table [X[ We 
use 'unsmoothed' values of ri/a which are the results of 
an independent fit to the heavy quark potential on each 
ensemble. Figure 18 shows the results for ri from each 
method as a function of the lattice spacing. Differences 
are evident on the very coarse lattices as a result of dis- 
cretisation errors but there is clear convergence as a ^ 0. 
The results are plotted against (a/ri)^ since the leading 
tree- level discretisation errors are at a^. Note that the 
behaviour of this plot is rather different from that ob- 
tained previously on the 2+1 flavor configurations (Fig- 
ure 3 of [21 ). There is a little less variation with a, to be 
expected because of the various improvements to the dis- 
cretisation of QCD. The main difference however is the 
direction of approach to a = 0. The results for ri from 
fr^^ are now very flat and the results from the T approach 
a = from below. This reflects a change in the relative 
discretisation errors of the quantities involved. 

As discussed above, the chiral fits involving the rjs give 
a result for ri of 0.3209(29) fm. We can fit the results 
from using the T method to test if they are consistent 
with this. Using the r]s result as a prior for the T fit 
then enables us to extract an improved result for ri which 
combines both methods. 



allowance for variations as a function of the sea light 
quark masses, although it is clear from the results that 
any such dependence is very small. Indeed the quanti- 
ties being used were chosen for their insensitivity to such 
effects. We also include an allowance for discretisation er- 
rors, both of the standard type, varying as {aA)^ and the 
from various am5-dependent type coming from radiative 
corrections in NRQCD. We also allow for the NRQCD 



systematic error from Tables VI and VII as a correlated 
error for all 5 ensembles. 

The fit to the T values using a large width prior for 
the physical value (0.32(10)fm) gives xV^of = 0.79 for 
5 degrees of freedom and ri(T) = 0.310(6) fm. This 
shows the required consistency in the determination of 
the lattice spacing from the two methods as a ^ 0. The 
fit including the prior value from the r]s analysis gives 
X^/dof = 0.76 and result: 



ri = 0.3209(26)fm. 



(32) 



This is slightly improved over the 77^ value on its own. 

This final value for ri can now be used to determine 
a on other ensembles if values of ri/a are available. We 
include in Table |X] the lattice spacing values on sets 1 to 
5 from using ri. 

Our result for ri can be compared to our previous re- 
sult of 0.3133(23) fm on the MILC 2+1 flavor dynamical 
asqtad lattices [21 . This is 2% lower than our current 
result with a combined uncertainty of 1% and so is not 
significant. In principle the two results do not have to 
agree because we are now including c quarks in the sea. 
However we expect this to have a small effect and then 
only in short-distance quantities [14 . We can obtain es- 
timates for the effect on the T 2S — IS splitting from 
the fact that it is proportional to the hyperfine split- 
ting. Missing c quarks in the sea increases the T mass 
by approximately 5 MeV with a smaller amount for ex- 
cited states. It therefore reduces the 2S — IS splitting 
by approximately 2.5 MeV or 0.4%. This could have 
led previously to a 0.4% underestimate of ri from the 
T 2S — IS splitting if ri/a itself was not affected. This 
effect is no larger than other sources of systematic error 
in the earlier calculation coming from radiative cor- 
rections to terms that are also now included. Thus we 
cannot claim to see any strong evidence of an effect from 
c quarks in the sea. Indeed if we compare the ri values 
coming from the T analysis alone there is a change of 
2(2)%, in which a 0.4% effect from c in the sea would be 
invisible. Any allowance for an effect on ri itself, also a 
fairly short distance quantity, would reduce this expected 
variation further. The rjs analysis would be expected to 
be very insensitive to sea charm because of the low in- 
ternal momenta inside these light hadrons. For that case 
we see only a 0.5% change in the value of ri obtained, 
again with a 2% error. 
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renormalisation. This latter renormalisation is given by: 



0.025 



FIG. 19: Values for the ratio of the b quark mass to the 
s quark mass in the MS scheme at a given scale plotted 
against the square of the lattice spacing. Results are ob- 
tained from combining NRQCD b quark masses and HISQ 
s quark masses with an 0{as) perturbative renormalisation. 
The errors on the points include statistical/fitting errors, lat- 
tice spacing errors and NRQCD systematic errors. The final 
result, including the 0(al) perturbative error is plotted as the 
shaded blue band. The result from our previous fully nonper- 
turbative calculation on ensembles including 2+1 flavors of 
sea quarks [481 149^ is given by the black filled circle at a = 0. 



VI. mb/rris 

From Table Ixnl we can determine the ratio of the bare 
NRQCD b quark mass to the bare HISQ s quark mass on 
each ensemble. To do this we must use the same determi- 
nation of the lattice spacing for the tuning of each mass, 
and so we use the lattice spacing determined from the T 
25'— IS splitting (columns 2 and 3). The lattice spacing 
error appears doubled in and once in 7715 because of 
their different dependence on the meson masses used to 
fix them. These errors are correlated in the ratio m^/m^ 
so one factor of the lattice spacing error cancels between 
numerator and denominator. 

The ratio of masses in different schemes (NRQCD and 
HISQ) is not particularly useful. However, we can con- 
vert this using perturbation theory to a ratio of masses in 
the same mass-independent scheme, such as MS, at the 
same scale, /i. The ratio then becomes scale-independent 
and the same in any scheme related to MS by a simple 
renormalisation. For both the NRQCD and the HISQ 
actions the mass renormalisation is known to ©(a^). 

The lattice to MS mass renormalisation constant is 
calculated by multiplying the lattice bare mass to pole 
mass renormalisation by the continuum pole mass to MS 



^(m) ^q,pole ( 1 + 



q,pole 



(33) 

The lattice bare mass to pole mass renormalisation for 
HISQ quarks is given for small quark masses by [50l 
M\ EH: 

ms,poie = ^ ( l + a,[- -In am, + 0.5387]... ) , (34) 



a 



TT 



where we have written the equation explicity for the 
strange quark mass. When equations [33] and [34l are com- 



bined to obtain the conversion factor from the lattice bare 
mass to the MS mass at scale ja and 0{as) the logarithm 
multiplying ag becomes ln(a/i), and there is a constant 
given by 0.5387-4/(37r). 

We can also write the NRQCD mass renormalisation 
in the form 



^b,pole 



ami) 



l + a,[--lnamb + .l^^«^^]... 



(35) 

although no In(am) term is explicit in that calculation. 
On doing this we find that the remainder term, ^Ni?QCD 
given in Table pCXIII[ has very little am^ dependence. 

Combining equations 33j 34 and 35 it is then clear that 
the ratio of MS masses for b and s is given to 0{as) by: 



-.MS 



(m) 



m 



MS 



ami) 
am. 



(36) 

where the ji dependence cancels out. The ratio of bare 
lattice masses from columns 2 and 3 of Table IXIII varies 
very little with lattice spacing with values between 51 and 
52. The renormalisation in equation [36] is a relatively 
mild one, with ag coefficient varying between 0.31 and 
0.39 with ami) value. We apply this one- loop renor maliza- 
tion with values taken as ay(1.8/a) from Table XXII 
The energy scale for as is then in agreement with the 
Brodsky-Lepage-Mackenzie scale calculated for the light 
quark (asqtad) mass renormalisation in ^Oj- This gives 
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the values for the MS mi) / ratio plotted in Figure 

The results in Figure [To] show very little dependence on 
lattice spacing or sea quark mass within the 1% statis- 
tical and systematic errors from the lattice calculation. 
A much larger error is that from missing higher order 
powers of ag in equation 36 We take account of this 



error by allowing a correlated error between the points 
of 1 X ay(1.8/a)^ along with a possible variation with 



ami) of the form ay(1.8/a)^ x (5x^/4 (see equation [17 
for a definition of 5xm)' This allows the term to have 
both a coefficient and a mass dependence which is three 
times that of the known term. We allow for possible 
dependence on sea quark masses and the lattice spacing 
by using a fit of the same form as that in equation 17 
The final fit result is then: 



-.MS 



m 



MS 



= 54.7(2.5), 



(37) 
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plotted as the shaded blue band in Figure [T9j The er- 
ror is dominated, not surprisingly, by the error from the 
unknown term. 

We can compare this new result to a combination of our 
earlier results for mij/mc (4.51(4)) from [48^ and mc/rris 
(11.85(16)) from [49 . These results were obtained en- 
tirely nonperturbatively by using the HISQ action for all 
the quarks. Then the ratio of lattice bare quark masses 
in the continuum limit is the ratio of MS masses at 
a given scale - the renormalisation factor cancels com- 
pletely. From the numbers above we have m^/mg = 
53.4(9) which is plotted as the black point at a = 
on Figure Our new, completely independent, result 



agrees well with this earlier value although it is much less 
accurate. 



VII. CONCLUSIONS 

In this paper we have determined the T spectrum using 
the NRQCD formalism for the b quarks in lattice QCD. 
We include several improvements over our earlier work. 
The key improvements are: 

• we use gluon field configurations with a fully 
0{asa'^) improved gluon action and HISQ quarks 
in the sea, provided by the MILC collaboration; 

• c quarks are now included in the sea; 

• we take the NRQCD action to a new level of accu- 
racy by including radiative corrections to the terms 
at next-to- leading relativistic order {v"^); 

• we improve the method for tuning the b quark mass 
so that systematic errors are reduced to 0.5%. 

With significantly improved systematic errors from 
NRQCD we are then able to determine the lattice spacing 
to better than 1% from the 25* — IS splitting. Using this 
we obtain M{hb) - M{1S) to 1.4% and M{V') - M(T) 
to 2.4% which is a strong test of NRQCD. This gives 
M{hb) = 9905(7) MeV to be compared with the exper- 
imental result of 9898.3(1.5) MeV [38] and M{V') = 
10375(22) MeV to be compared to the experimental re- 
sult of 10355.2(5) MeV [35]. 

We have examined the T and 775 dispersion relations 
in much more detail than before, so that we can quantify 
the effect of the radiative corrections to the kinetic 
terms in the action. We are also able to show how small 
are the deviations from continuum rotational invariance. 
This enables us to tune the b quark mass to 0.5%. 

Our result for the hyperfine splitting between the T 
and 775 states is much more accurate than in our earlier 
work because we have included the critical renormalisa- 
tion of C4 (the coefficient of the a-B term ) in our analysis. 
We obtain M(T) - M{r]b) = 70(9) MeV now with a 13% 
error. This gives M{r]b) of 9390(9) MeV to be compared 
to the experimental result of 9390.9(2.8) MeV. 
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FIG. 20: The spectrum of bottomonium states from lat- 
tice NRQCD (colored symbols with error bars) compared to 
experiment (black lines). Blue crosses denote results used to 
tune parameters, pink open squares results to be compared to 
experiment and red circles predictions ahead of experiment j^] 

^ For simplicity we mark the T with a blue cross although in fact 
we use the spin-average of T and 775 to tune the b quark mass. 



Our result for M(T') — M{iri'^) is also much more ac- 
curate largely because of a huge improvement in the sta- 
tistical error. We find a 2S hyperfine splitting that is 
half as big as the IS hyperfine splitting at 35(3) MeV, 
predicting M{r]'^) = 9988(3) MeV. 

These new results are collected together in a plot of the 



T spectrum from improved lattice NRQCD in Figure 20 



We mark with different symbols those results used to 
tune parameters, those which correspond to masses al- 
ready known from experiment, and those (the r]'^) which 
are predictions. We include the P-wave fine structure 
from our results for C4 = 1.15 on the fine lattices, set 5, 
since this C4 is close to the perturbative value on those 
lattices. We include an additional 10% error for missing 

terms in our NRQCD action. D-wave T masses from 
our calculation will be reported elsewhere. 

Light meson (tt, K and rjs) masses and decay constants 
are also given here that enable us to determine the prop- 
erties of the rjs nieson and give a complementary deter- 
mination of the lattice spacing to better than 1%. The 
calculation shows significantly improved discretisation er- 
rors over our earlier results on ensembles including 2+1 
flavors of asqtad quarks [21 . The results on the prop- 
erties of the are in agreement with our earlier work. 
However, our earlier result was not able to distinguish 
the mass and decay constant of the r]s from that that 
would be obtained in leading order chiral perturbation 
theory. We now obtain M{r]s) = 0.6893(12) GeV and 
/r^^ =0.1819(5) GeV. In both cases these values disagree 
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significantly, but by less than 1%, from the leading or- 
der expectation. The r]s particle is relatively insensitive 
to sea u/d quark masses and so it is very useful to have 
accurate results for its properties for tuning the s quark 
mass and determining the lattice spacing on other lattice 
ensembles. 

The T 25* — IS* and rjg determinations of the lattice 
spacing can be compared through a third parameter, ri, 
from the heavy quark potential. We show that both de- 
terminations agree in the continuum and chiral limits and 
give a physical value for ri of 0.3209(26) fm. This can 
also be used to determine the lattice spacing on other 
lattice ensembles. 

We also combine T and r]s calculations through a de- 
termination of the ratio of MS b quark to s quark masses 
of mb/rris = 54.7(2.5), in agreement with our earlier re- 
sult from HISQ quarks alone of 53.4(9) [48l[49]. 

Finally we comment on the effect of including c quarks 
in the sea. We have seen no significant effect on any of 
the observables that we have calculated compared to re- 
sults obtained with 2+1 flavors of sea quarks. The results 
cannot be compared lattice spacing by lattice spacing be- 
cause of changes to the lattice QCD action that reduce 
the size of discretisation errors in our new results. Fi- 
nal physical results can be compared, however, with and 
without sea c quarks to see if there is a difference. In our 
earlier 2+1 flavor calculations we estimated that the 
presence of sea c quarks would shift the T and 7^5 masses 
downwards by 5 MeV through an induced additional local 
potential which was proportional to a^(5^(r)/m^. This 
would have a smaller effect on radial excitations of the 
T than on the ground state masses and very little ef- 
fect on P-wave states. We then estimate the effect on, 
for example the T IP - 15 splitting to be (9(1%). This 
would barely be visible above the errors in our current 
calculation and the errors in the earlier calculation were 
somewhat larger, so any comparison certainly has an er- 
ror of greater than 1%. However, it is clear from our re- 
sults that no unexpectedly large effect has appeared. For 
light hadrons we expect even smaller effects and there we 
can limit any differences in M^^ and /^^ to smaller than 
1%, with the main error coming from our earlier calcula- 
tion ^21^. 

We are now combining b quark propagators from our 
improved NRQCD action with /, s and c propagators on 
these ensembles to study P, P5 and Be meson masses 
and matrix elements. Significantly improved systematic 
errors should be possible both from the NRQCD action 
and because we are working much closer to physical light 
sea quark masses than before with an improved gluon 
and sea quark action. 
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Appendix A: Gauge Action 

For clarity, the gauge action So used in the genera- 
tion of the MILC ensembles will be summarised in this 
section. See [19 . The action is a tadpole and one- loop 
improved Liischer-Weisz action. 



Sg = P 



ReTr(i?) 



R ^ 

+ CTE(l-^I^eTr(T)^ 



(Al) 



where the sums are over plaquettes P, rectangles R 
and twisted loops (or parallelograms) T. The coeffi- 
cients are calculated perturbatively through 0{as) in- 
cluding both gluonic loops [52 and contributions from 
HISQ sea quarks [20 . The tadpole improvement param- 
eter was chosen to be the fourth root of the plaquette 
uqp = (|Re Tr(P))^/^ and, via a perturbative calcula- 
tion of the plaquette, gives an expression for the strong 
coupling constant ag = —1.303615 log iiop- uqp also ap- 
pears in the gauge coupling as /3 = 10/{g'^UQp). The 
coefficients used are 



Cp = 1.0 

Cr 



(1 - (0.6264 - 1.1746A^/) log(iigp)) 



Ct = ^(0.0433 -0.0156A^/)log(^xgp) (A2) 

Uqp 

The inclusion of these terms mean that the gauge ac- 
tion is improved completely through order 0{asa'^). As 
mentioned in the text, sea quarks are included using the 
HISQ action [16] with a U{3) projection (only) for the 
intermediate re-unitarization step. 



Appendix B: Perturbative determination of 
radiative corrections to Ci coefficients in the 
NRQCD action and the mass renormalization 



We are grateful to the MILC collaboration for the use 
of their gauge configurations and particularly to Doug 
Toussaint for help in reading them and in providing val- 



Spin-independent coefficients. The q coefficients ap- 
pearing in the NRQCD action, equation [2j have ex- 
pansion 1 + + O(q^s). The ^ for the kinetic 
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FIG. 21: The Feynman diagrams needed for the calculation 
of 0{as) corrections to the heavy quark self energy. From top 
to bottom, rainbow, tadpole and uq counterterm diagrams. 



terms, i = 1,5,6, are determined following the method 
of [26l[53]. The NRQCD quark self-energy is calculated 

through 0{as) and the c-^^ are given by the require- 
ment that the correct energy-momentum relationship be 
obtained through 0(asV^). The terms proportional to 
(A'^^^)^ in equation |2] can be merged together so that 
this term in 6H appears as: 



(A(2))2 



am\) 



(Bl) 



Thus only two radiative corrections need to be calculated 
for the complete set of kinetic terms at 0{y^)^ i.e. for 
ci and C5. The radiative correction for ci then applies 
equally to Ci and cq. 

The full inverse NRQCD quark propagator at 0{as) is 



aG-^{p) = Q-^{p)-asaT.{p) 



(B2) 



where ap = {ap^ap^) is a 4- vector in lattice Euclidean 
space. The pole in the propagator is identified as aco'(p) = 
ip^a. The expansion of co'(p) in powers of the spatial mo- 
mentum can be used to identify the quark mass renor- 
malisation factor and wavefunction renormalisation 
factor Z2 but also to tune c^^'^ and €5^'^ to appropriate 
values. Q~^{p) is the quark propagator obtained at tree 
level from the NRQCD action, including the (as yet un- 
known) radiative corrections to ci and C5 . Its pole is then 
given by: 



acjo(p) 



2am5 



8{am}))^ 
1 1 



2n ami) 



r (1) a^p 
' \ ^ 24ar7 
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24am5 



(B3) 



^{p) is the one-loop self-energy and consists, as shown 
of rainbow and tadpole diagrams as well 



in Figure 21 



as diagrams containing insertions of the one-loop piece 
of the tadpole- improvement factor, uq. Writing uq = 



1 + asU^\ we have Uq^"^ = 0.750 for the Landau link 
tadpole parameter, uql [54 . We have 



cj(p) = a;o(p) - as^{uJo{p),p) 
and can expand S to as: 



a^{p) = SoH + Ei(cj) 



2am}) 



(B4) 



(B5) 



S2(w) 



The are extracted from suitable combinations of par- 
tial derivatives of S: 



So 
Si 



aS(p = 0) 

arrib- 



(B6) 
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Each of the also has an expansion in powers of uj as 
Yi 



S^oSf^cj^ Then 



aco'(p) 



a2p2 



P' 



2ami)^r S{ami)^r)^ 
where m^^r = Zmrrib and 

Z„ = l + a,Z« =l + a,(s(') 



asa5uj{p) (B7) 



to this order. The correction term Suj is given by: 
aSuj = Wo ^ 



(B8) 
(B9) 
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(BIO) 



The requirement that lattice NRQCD reproduce the low- 
energy physics of full QCD means that 5uo can only be a 
pure energy shift independent of spatial momentum, i.e. 
the coefficients of (p^)^ and p^ in equation BIO must be 
zero. Thus 
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(Bll) 
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The Feynman rules were generated automatically us- 
ing the HiPPy package and the Feynman lattice integrals 
for S and its derivatives were constructed and evaluated 
numerically using the HPsrc package and VEGAS con- 
tained therein [55l |56] . We use analytic differentiation 
using the TaylUR package in the HPsrc Fortran code 
together with numerical differentiation which, for suffi- 
ciently smooth functions, can be up to an order of mag- 
nitude faster than analytic differentiation. 

Because the kinetic (a^p^)^ term is included at tree 
level both Wi and W2 are Infra- Red (IR) finite and so no 
IR regulation is needed although a gluon mass was used 
to regularize intermediate divergences. However, the in- 
tegrals arising from the rainbow diagram still have large 
peaks in the IR region. These peaks arise because the 
differentiation generates extra powers of the heavy quark 
NRQCD propagator in the integrand. In this case to use 
numerical differentiation alone proves to be unstable and 
it is imperative to use a mixture of analytic and numerical 
approaches and also to introduce a suitable subtraction 
function to remove the most severe behaviour of the inte- 
grand. In contrast, the integrals arising from the tadpole 
diagram are well behaved because they contain no quark 
propagators but they are expensive to evaluate since the 
two-gluon vertex contains a large number of terms. In 
this case, numerical differentiation proved to be the most 
efficient for the higher order mixed derivatives without 
compromising accuracy. In all cases the temporal deriva- 
tives were done using the analytic method. 

We checked that the results agree well with those of 
Morningstar [26l |53] for his gluon and NRQCD actions. 
For the simplest gluon and NRQCD actions the results 
agree with the analytic calculation of Monahan |57] . 

The contribution to the c,-^^ from the u^"^ insertions of 



Figure 21 can be calculated analytically. This gives: 
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These contributions are sizeable and act to cancel contri- 
butions coming from the other diagrams, as part of the 
'tadpole- improvement' mechanism [5j [26]. This is par- 



ticularly true for c^^-'; less so for c^^^, as in [26|. The 

cl^^ values will then change depending on the tadpole- 
improvement parameter chosen, for example uqp or 
uoL, because the q must compensate perturbatively for 
changes in uq. Here we use uql in the NRQCD action 
and this is the only uq that affects the q to 0{as). uqp 
is used in the gluon action and counterterms from this 
will app ear in the q at higher order. 

Table XXI gives the results for Ci'^ and c^^ for 3 
different values of am^ and stability parameter n = 4. 
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TABLE XXI: Coefficients c[^^ and that multiply as in 
the one-loop correction to the kinetic terms in the NRQCD 
action used here in conjunction with the improved gluon ac- 
tion described in Appendix [A] 



arrib 


n 




^5 


1.95 

2.8 

3.4 


4 
4 
4 


0.774(21) 
0.951(26) 
0.952(30) 


0.392(17) 
0.406(11) 
0.445(10) 


TABLE XXII: Values for ay used in calculating 
corrected coefficients ci, C5, cq and C4. 


the 0{as)- 


Sets 


1/a 
GeV 






fine 
coarse 
very coarse 


2.2 
1.6 
1.3 


0.32 0.28 
0.39 0.33 
0.46 0.38 


0.225 
0.255 
0.275 



The mass values are not exactly the am^ values used for 
the numerical work in this paper but the c-^'' show very 
mild dependence on am^, so we can simply interpolate 
to the ami) values we are using. The results are differ- 
ent from those of [26 because both the gluon action and 
the NRQCD action have changed. However, qualitative 
features are the same in that the values are not large 
and only mild dependence on am^ is seen for am^ larger 
than 1. We note that the c-^^ coefficients will change 
if higher order terms are added to the NRQCD action. 
For example [27 tests were done with an NRQCD action 
which included a term in SH of —A'^^^/ (ISOam^), remov- 
ing 0{a^) discretisation errors from Hq. This changes 
Cg^^ to 0.017(4) for am^ = 1.95 and n = 4 to be com- 



XXI 



pared with the result of 0.392(17) in Table 

The coefficients in Table [XXll need to be combined with 
a value for ag to give final results for ci, C5 and cq. The 
scale, g*, used for as was taken as that calculated for the 
Brodsky-Lepage-Mackenzie scheme in Figure 10 of [26] . 
assuming that this does not change significantly with the 
changes in the action used. This gives ~ 1.4/a for C5 
and g'* ^ 1.8/a for ci. We take as from [48 , specifically 
the value aj^{Mz^nf = 5) = 0.1183. We convert this 
to ay [53 and run perturbatively to values using Uf = A 
and appropriate scales q* . The q* values are calculated 
for the very coarse, coarse and fine ensembles using ~ 
1.3, 1.6 and 2.2 GeV respectively from Table [Xj The ay 
values obtained are listed in Table IXXHl These are com- 
bined with the coefficients in Table jXXl to give the values 
used in Table [ll[ The coefficient was reduced slightly 
on the fine lattices (to 0.766) to account for the fact that 
the value of am^ used was slightly smaller than that for 
which the coefficient was calculated. 

The remaining errors in the kinetic q coefficients after 
this one-loo p corr ection has been made will be 0{a^g). 
From Table XXII we can see that O.Sofy ranges from 0.1 



on the very coarse ensembles to 0.05 on fine. The impact 
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TABLE XXIII: Coefficients and 



that multiply as 
in the one-loop correction to the mass renormalization and 
the a • B term in the NRQCD action respectively. These 
are calculated with the NRQCD action used here and the 
improved gluon action described in Appendix |a1 ^nrqcd 

m 

column 4 is Zm 

+ 2 \n(amb) / TT , as described in the text. 



arub 


n 




^NRQCD 




1.9 


4 


0.439(3) 


0.848(3) 


0.691(7) 


2.65 


4 


0.263(5) 


0.883(5) 


0.775(8) 


3.4 


4 


0.150(3) 


0.929(3) 


0.818(4) 




FIG. 22: The 0{as) coefficient in the perturbative expansion 
of C4, coefficient of the a • B term, plotted against the bare b 
quark mass. 



of these errors can then be estimated from the size of the 
effect that we see on physical observables from the 0{as) 
corrections. 

We have not corrected the coefficient of the Darwin 
term, C2, in our NRQCD Hamiltonian. We have however 
assessed the effect of taking C2 to be as large as 1.25 on 
the meson kinetic mass and on the hyperfine splitting 
and we find the effects to be small. 

Mass renormalization. In Table [XXTTT] we give values 

in the mass renormalization 
Zjn was calculated previously for differ- 
ent parameter values in [59 . We also show the result 

(i) ^NRQCD ^NRQCD ^^^i 



for zin^ , the coefficient of a 



of equation B8 



of adding 2 ln(am5 )/7r to Zrr 
be used in section VI to derive the ratio mi^/ms in the 



MS scheme. Only mild dependence on amj) is seen in 

^NRQCD 

Spin- dependent terms. The perturbative renormalisa- 
tion of field-dependent terms has to be done in a different 
way and this includes all the spin-dependent terms. Re- 
cently the radiative corrections to C4 have become avail- 
able [28] . They were calculated by matching the effective 



TABLE XXIV: Coefficients di and d2 of spin-dependent 4- 
quark operators that give rise to a correction to the hyperfine 
splitting. These are calculated with the NRQCD action used 
here with the parameters given in columns 1 and 2 and the 
improved gluon action described in Appendix [A] 



amb 


n 


di 


d2 


1.9 


4 


-ln(1.9) + 0.796(4) 


ln(1.9)/3 -0.311(1) 


2.65 


4 


-ln(2.65) + 0.448(6) 


ln(2.65)/3 -0.195(2) 


3.4 


4 


-ln(3.4) + 0.038(8) 


ln(3.4)/3 -0.058(2) 



action in NRQCD to continuum QCD using the back- 
ground field method. We give values for €4^^ in Ta- 



^4 



ble |XXIII| appropriate to the am^ and n values we are 
using here. 

A number of pieces go in to the calculation of 
These include renormalisation of the chromomagnetic 
moment, renormalisation of the wavefunction and, be- 
cause C4 multiplies a term in the bare lattice NRQCD 
Hamiltonian which includes the bare lattice quark mass. 



the mass renormalisation. 



is the sum of two 



pieces, one of which has polynomial dependence on 
the bare quark mass amb and the other is propor- 
tional to log (amfo). The logarithmic term has coefficient 
—3/ (27r) 28|. Both terms are included in the total result 
given in Table pCXIII| - the logarithmic term is of similar 
size to the polynomial term over the range of am^ that 
we are using. The c^^-' valu es are combined with values 



for ay\7r/a) (given in Table XXII) to give the results for 



the 1-loop corrected C4 coefficients given in Table pOTT 



Figure 
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gives a more complete picture of c^^^ by plot- 



ting values as a function of am^. We see relatively little 
ami) dependence until am^ becomes smaller than 1 when 
€4^^ starts to diverge. This is typical of the behaviour 
of radiative corrections to coefficients in the NRQCD ac- 
tion. 

Finally we give results for the coefficients of the spin- 
dependent 4-quark operators that contribute to the hy- 
perfine splitting. These terms have coefficients di and ^2 
that multiply terms that would appear in the NRQCD 
action at order 0{alv^): 



di 



{am^Y 



' [ami))- 



r(V^W)-(x^^V^). (B13) 



These terms are subleading compared to treelevel op- 
erators but contribute to the hyperfine splitting at the 
same order as as corrections to C4 [28 . We do not in- 
clude these 4-quark terms in our NRQCD action but we 
can estimate their effect because they give rise a shift in 
the relative energies of the T and 775 which is proportional 
to the 'wavefunction-at-the-origin', and given by: 



6al{d2-di) 



mf 



(B14) 
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The relevant coefficients, di and d2 are given in Ta- 
ble |XXIV[ They were calculated previously for a slightly 
different NRQCD action in [28]. The coefficients are re- 
lated by: 



di = -3d2- -{1-H2)) 



(B15) 



where the term proportional to (In 2 — 1) is from 775 an- 
nihilation to 2 gluons. This term increases the hyperfine 
splitting from equation B14| since it pushes the 7^5 mass 



down. When equation B14 is applied for just this piece 
we obtain an estimate of its size of about 1 MeV. This 
is smaller, but in agreement with, the earlier estimate 
of 2.4(2.4) MeV applied in section IIIC in deriving an 
appropriate spin- averaged l^* meson mass to tune the b 
quark mass against. For that purpose the shift is com- 
pletely negligible, representing a tiny fraction of the 775 
mass. For the hyperfine splitting it is a more important 
issue. For this we take the results from equation |B14 
because that provides a consistent treatment of all the 
4-quark operator effects. 

di and ^2 are separated into logarithmic and nonloga- 
rithmic pieces in Table pDCIV The nonlogarithmic piece 
has significant mass dependence here, becoming small at 
large ma. The logarithmic and nonlogarithmic terms in 
fact cancel for ma around 1.9. In assessing the error in 
the hyperfine splitting from missing higher order terms 
multiplying the 4-quark operators we are careful not to 
assume that this is generic behaviour. 

We combine th e di an d ^2 coefficients with ay(7r/a) 
values from Table XXII and values for |?/;(0)p from our 
fits to obtain corrections to the hyperfine splitting that 



are applied in section HIE 3 



Appendix C: Nonperturbative determination of 
radiative corrections to C3 and C4 coefficients in the 
NRQCD action 

An alternative approach to tuning the spin-dependent 
coefficients C3 and C4 is from matching fine structure 
in the spectrum to experiment. Here we use the P- 
wave fine structure in the T spectrum to do this O [12] . 
The P-wave masses are shifted from the spin-average by 
amounts that depend on spin-spin coupling terms propor- 
tional to S • S or Sij and spin-orbit terms proportional 
to L • S. The spin-spin terms are proportional to C4 and 
the spin-orbit terms to C3. We can make a combinations 
of the state masses in which the eigenvalues of Sij 
and S • S cancel, and a separate combination in which 
the eigenvalues of S • S and L • S cancel. Thus each of 
these combinations gives a result which should depend 
on only one of the spin-dependent couplings in our cur- 
rent NRQCD action. Comparing these combinations to 
experiment allows us to tune C3 and C4. Note that 4- 
quark operators, discussed in Appendix [B] with reference 
to their effect on S'-wave hyperfine splittings, have neg- 
ligible effect on P-wave states, and they therefore give a 
very clean determination of C3 and C4. 
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FIG. 23: The masses of the lowest-lying P-wave states in the 
T spectrum plotted relative to the spin- average of the ^P 
states for the coarse lattices, sets 3 and 4 (top plot) and the 
fine lattices, set 5 (lower plot). In each plot we compare mass 
splittings for C4 = 1 with a nonperturbatively tuned value 
for C4 chosen to match a combination of mass splittings to 
experiment (see text). C3 = 1 in all cases. For the ^P2 states 
the ^Pe is plotted to the left of the ^Pt2- 



The eigenvalues for L-S, Sij and S-S for '^^'^^Lj states 
{3Po,^Pi,^P2,^Pi} (i.e. {x60,X6i,X62,/^6}) are: 



L S 

Sij 
S-S 



{-2,-1,1,0}; 

{-4,2,-2/5,0}; 

{1/4,1/4,1/4,-3/4}. 



(CI) 



We see that the S • S terms affect the splitting between 
the ^Pi and the spin-average of the ^P states. We ex- 
pect this splitting to be small because, in a potential 
model approach, the accompanying spin-dependent po- 
tential would be a delta function at the origin with very 
little overlap for P-wave states. Both Sij and L • S terms 
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aE 
aE 
aE 
aE 
aE 



aE 
aE 
aE 
aE 
aE 



I'Pi) 
l^Po) - 
l^Pi) - 
I^Pe) 
I'Pt.) 



aE{^P) 
aEC'P) 
aE[^P_) 
- aE(^P) 
- aEC'P) 



Set 3 

arub = 2.66 
C4 1.0 
C3 = 1.0 



Set 3 

arub = 2.66 
C4 1.25 
C3 = 1.0 



Set 4 

amb = 2.62 
C4 = 1.25 
C3 = 1.0 



Set 5 

amb = 1.91 
C4 = 1.0 
C3 = 1.0 



Set 5 

amb = 1.91 
C4 = 1.15 
C3 0.96 



0.5655(23) 
0.5460(20) 
0.5611(24) 
0.5747(30) 
0.5732(28) 



0.5247(22) 
0.5017(20) 
0.5213(26) 
0.5359(28) 
0.5331(28) 



0.5253(20) 
0.5034(18) 
0.5218(24) 
0.5354(25) 
0.5328(27) 



0.4833(10) 

0.4678(9) 

0.4802(10) 

0.4903(11) 

0.4893(11) 



0.4478(11) 

0.4312(9) 

0.4454(11) 

0.4549(12) 

0.4538(12) 



-0.0010(7) 

-0.0204(13) 

-0.0053(8) 

0.0082(12) 

0.0067(8) 



-0.0016(10) 

-0.0246(14) 

-0.0050(10) 

0.0096(12) 

0.0068(10) 



-0.0011(9) 

-0.0231(13) 

-0.0046(9) 

0.0089(11) 

0.0064(10) 



-0.0008(2) 

-0.0163(4) 

-0.0039(2) 

0.0062(3) 

0.0052(3) 



-0.0009(2) 

-0.0176(5) 

-0.0033(3) 

0.0062(4) 

0.0051(3) 



2aE{i^P2E) + 3aE(i^P2T2) 


0.093(8) 


0.104(9) 


0.097(8) 


0.072(3) 


0.073(3) 


-3aE{l^Pi) - 2aE{l^Po) 












0AaEll^P2E) + 0.6aE{l^P2T2) 


-0.018(2) 


-0.026(4) 


-0.025(3) 


-0.0153(7) 


-0.0198(9) 


-3a^(l^Pi) + 2a^(l^Po) 













TABLE XXV: Fitted energies for P-wave states on sets 3, 4 and 5 for the NRQCD parameters given. The last two rows give 
the combination of energies used to determine cs and C4. Errors are from statistics/fitting only. 



3 3 4 5 5 

C4 = 1 C4 = 1.25 C4 = 1.25 C4 = 1 C4 = 1.15 

C3 1.0 C3 1.0 C3 1.0 C3 1.0 C3 0.96 



5^(1^P2) 
^(l'P2) - 



- 3£^(rP2) - 2^(rP2) 0.151(13) 0.168(14) 0.160(13) 0.161(6) 0.162(7) 
3^(1^P2) + 2^(1^P2) -0.028(4) -0.042(6) -0.041(5) -0.0342(16) -0.0442(20) 



TABLE XXVL Combinations of P-wave energies needed to fix cs and C4 in GeV. Lattice spacing values are taken from Table pC] 
Errors are statistical/fitting only. 



affect the splittings within the ^P sector but not the split- 
ting between ^Pi and ^P. The L • S terms give the con- 
ventional ordering of ^Pq, ^Pi and ^P2, but with the ^P2 
splitting from the ^Pi larger than that between the ^Pi 
and ^Pq. The Sij terms will push down the ^P2 relative 
to the others. Thus the final sphttings depend on the rel- 
ative strength of the accompanying potentials for these 
terms and, in NRQCD language, the coefficients C3 and 
C4. 

The combination of spin-splittings that depends on C4 
(through Sij) and is independent of C3 is |2] 



M(x62)-3M(x6i) + 2M(x6o) 



(C2) 



with experimental value: -47.4(1.3) MeV [35 , determin- 
ing the errors on mass differences by adding the errors 
on the masses in quadrature. Likewise the combination 
that depends on cs only 



5M(X62) - 3M(x6i) - 2M(x6o) 



(C3) 



with experimental value: 163.8(3.2) MeV [35 . Compari- 
son of our results with experiment for these combinations 
then allows us to fix cs and C4. 



Table XXV gives the results for the energies in lattice 
units of the lowest-lying P-wave states for the amb val- 
ues given in Table III on the coarse (sets 3 and 4) and 



exponential fits of the form given in equation [7] to the 
2x2 matrix of correlators for each P-wave meson done 
as a single simultaneous fit. This enables us to extract 
mass differences more precisely from the fit than the in- 
dividual masses and the splitting between each state and 
the spin-average of all the ^P states is also given. Note 
that we give separate values for ^P2E and ^P2T2 lattice 
representations of the J = 2 state. Differences between 
the values obtained for T2 and E would be a sign of dis- 
cretisation errors. We do not have a significant signal for 
this but the E state is higher than the T2 in all of our fits. 
The difference is about 3(2) MeV on the coarse lattices 
and 2(1) MeV on the fine lattices. 

Results are given for the case C4 = 1 and a nonzero 
value of C4 chosen to give reasonable agreement with the 



experimental value of the combination in equation C2 



Results in GeV for the two combinations tested are given 
in Table pCX VI using values of the lattice spacing from the 
(25* — IS*) splitting in Table |xj and combining the E and 
T2 representations for the ^P2 state with the appropriate 
number of spin states. We see that C3 = 1 within our 
errors, but C4 needs to be larger than 1, more so on the 
coarse lattices than the fine. We take the same value of C4 
on both coarse sets since the tuning should not depend on 
the sea quark masses and indeed our results demonstrate 
that it does not. 



fine (set 5) lattices. The results were obtained from 5 



Figure 23 shows the spectrum of P-wave states relative 



36 



to the spin-average (5M(x62) + 3M(x6i) + M{xbo))- 
The results with C4 > 1 clearly agree better with experi- 
ment than for C4 = 1. The main effect of increasing C4 is 
to push the Xbo state down relative to the spin-average. 
Very little else changes. In particular we see that the 
splitting between the ^Pi and the ^P spin- average is very 
small and negative in all cases. It increases with increas- 
ing C4 but, because the splitting itself is so small, this is 
not significant. We obtain a P-wave hyperfine splitting 
of 2(2) MeV on both coarse and fine lattices, where the 
experimental result is 1.6 ± 1.5 MeV [38] . 

On the fine lattices we took C4 = 1.15, based on ini- 



tial calculations. From Table XXVI and Figure 23 this 
appears to be an underestimate and an improved value 
would be 1.18. We also used C3 = 0.96 but that value 
is indistinguishable within our errors from C3 = 1.0. On 
the coarse lattices, C4 = 1.25 is also a slight underesti- 
mate, although it agrees within statistical errors with the 
correct answer. 

The final nonperturbatively tuned values for cs and C4 
that we obtain are then: C3 = 1.00(9) (2) (10) on coarse 
lattices and 1.00(4) (2) (10) on fine lattices. Our best es- 
timates for C4 are: 



C4(coarse) = 1.28(7)(1)(5) 
C4(fine) = 1.18(2)(1)(5). 



(C4) 



The first error is from the statistical/fitting error on the 
P-wave masses along with the lattice spacing error. The 
second error is from experiment. The third error is a sys- 
tematic error from terms in NRQCD that are missing 
from our calculation but that have effectively been ab- 
sorbed into the value of C3/C4 from matching to experi- 
ment. The spin splittings could change by (9(10%) from 

terms. Note that C4 gets closer to 1 as the lattice gets 
finer as we expect. Note also that the relationship be- 
tween Cs and C4 that holds in potential models or NRQCD 
in the continuum because of Lorentz covariance [36^ is not 
applicable to lattice NRQCD in this formulation. 

The agreement between the C4 coefficients obtained 



nonperturbatively and the C4 coefficients obtained at 
0{as) in Appendix [B] is good, and certainly within pos- 
sib le v ariation of the perturbative coefficients (see Ta- 
ble 



XIII) 



Our nonperturbative results for C4 and C3 agree well 
with those derived for the same NRQCD action on differ- 
ent gluon configurations at similar lattice spacing in [12] . 
There spin-dependent terms at are also included in a 
separate calculation and then the values derived for the 
v"^ coefficients C3 and C4 change (both increasing). Note 
that in that paper the calculations were done with tree- 
level coefficients for all q and the results rescaled for the 
derived values of C4. 



Appendix D: Results for the Kinetic mass 

We give here the tables of results for the ground-state 
energies and kinetic masses (as defined in equation 10 ) of 



the T and 775 mesons for different lattice meson momenta 
in units of 27ra/L. We also give the spin-averaged IS* 
kinetic mass. Results are taken from simultaneous fits 
to local correlators with the given momentum and with 
momentum zero using the fit form in equation [7] The 
energies for each momentum and the energy difference 
which yields the kinetic mass are given directly by the 
fit. Correlations between the correlators mean that the 
error on the energy difference is typically smaller than the 
combined errors from the separate energies. This is par- 
ticularly true of the 'on-axis' momenta which have only 
one non-zero component. So, for example, the kinetic 
mass for momentum (2,0,0) is more precise than that for 

(1,1,1). 

The results are given for coarse ensemble set 3 and fine 
ensemble set 5 with separate results for the case where 
Ci,5,6 are taken to be 1 and the case where Ci^^^q are a^- 
improved. We take 9 exponential fits on set 3 and 7 
exponential fits on set 5. 
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(0,0,0) (1,0,0) (1,1,1) (2,0,0) (2,2,1) (3,0,0) 

aECSo,P) 0.25529(4) 0.26119(4) 0.27309(4) 0.27890(4) 0.30830(8) 0.30814(6) 

aECSi.P) 0.28626(6) 0.29220(7) 0.30426(7) 0.31007(9) 0.33977(17) 0.33957(14) 

aMKin(ryb) - 5.773(10) 5.767(7) 5.788(2) 5.787(7) 5.805(3) 

aMKin(T) - 5.716(25) 5.703(17) 5.739(7) 5.732(15) 5.753(11) 

aMKin(16') - 5.730(20) 5.719(14) 5.751(6) 5.746(12) 5.766(8) 



TABLE XXVII: T and r/b energies and kinetic masses in lattice units for various lattice momenta on coarse set 3 for b quark 
mass arrib = 2.66 and ci,5,6 set to 1. 



(0,0,0) (1,0,0) (1,1,0) (1,1,1) (2,0,0) (2,1,1) (2,2,1) (3,0,0) 

aE{'So,P) 0.26096(4) 0.26684(4) 0.27273(4) 0.27860(4) 0.28438(4) 0.29610(4) 0.31348(6) 0.31335(6) 

aE{^Si,P) 0.29243(6) 0.29838(6) 0.30434(6) 0.31030(7) 0.31611(8) 0.32799(8) 0.34555(14) 0.34536(14) 

aMKinivb) - 5.818(7) 5.819(7) 5.817(7) 5.839(3) 5.834(4) 5.844(7) 5.859(4) 

aMKin(T) - 5.747(18) 5.748(17) 5.742(17) 5.778(7) 5.764(9) 5.778(13) 5.798(12) 

^Kin(16') - 5.764(15) 5.766(14) 5.761(14) 5.793(5) 5.782(8) 5.795(11) 5.813(10) 

TABLE XXVIII: T and r/b energies and kinetic masses in lattice units for various lattice momenta on coarse set 3 for b quark 



mass arrib = 2.66 and ci,5,6 set to their 0{as) improved values. Slight differences with Table VIII for zero momentum energies 
arise because we fit a single zero momentum correlator rather than a 5 x 5 matrix. 



(0,0,0) (1,0,0) (1,1,1) (2,0,0) (2,2,1) (3,0,0) 

aE{'So,P) 0.24652(3) 0.25107(3) 0.26010(3) 0.26461(3) 0.28713(4) 0.28712(4) 

aE{^Si,P) 0.27153(5) 0.27610(4) 0.28518(5) 0.28974(5) 0.31244(6) 0.31246(7) 

aMKinivb) - 4.244(6) 4.252(6) 4.2548(14) 4.2516(23) 4.251(3) 

aMKin(T) - 4.222(14) 4.230(13) 4.225(3) 4.223(4) 4.215(6) 

^Kin(16') - 4.228(12) 4.236(11) 4.2327(24) 4.230(4) 4.224(5) 



TABLE XXIX: T and r/b energies and kinetic masses in lattice units for various lattice momenta on fine set 5 for b quark mass 
arub — 1.91 and ci,5,6 set to 1 



(0,0,0) (1,0,0) (1,1,0) (1,1,1) (2,0,0) (2,2,1) (3,0,0) 

aE{'So,P) 0.25827(3) 0.26278(3) 0.26727(4) 0.27173(3) 0.27620(3) 0.29850(4) 0.29847(3) 

aEl^Si.P) 0.28390(5) 0.28844(4) 0.29299(6) 0.29747(5) 0.30199(5) 0.32451(6) 0.32447(6) 

aM^J(r]b) - 4.278(7) 4.286(9) 4.287(6) 4.2914(14) 4.2920(23) 4.2951(19) 

aMKin(T) - 4.245(15) 4.251(17) 4.256(14) 4.2515(29) 4.2523(44) 4.2538(41) 

^Kin(16') - 4.253(12) 4.260(15) 4.264(11) 4.2615(24) 4.2622(37) 4.2641(35) 



TABLE XXX: T and rjb energies and kinetic masses in lattice units for various lattice momenta on fine set 5 for b quark mass 
arrib — 1-91 and ci,5,6 set to their 0{as) improved values. Slight differences in zero momentum energies are seen compared to 
Table [Vm] because we used i^ol =0.85246 rather than 0.8525 and are fitting to single correlators rather than a 5 x 5 matrix. 



